{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Simulation of a Kerr-Black-Hole\n",
"by **Florian Hollants**"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This Jupyter Notebook simulates a rotating Black Hole in SageMath using the Kerr Metric. The Project is based on the work of Florentin Jaffredo. The parameters for mass and angular momentum can be altered in Cell 5 by changing \"m_val\" and \"a_val\", while the \"Angle\" can be changed in Cell 18. "
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"# n_cpu = 4 # 4 Go Ram minimum\n",
"# n_geod = 100\n",
"# nx, ny = 180, 90"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"n_cpu = 8 # 8 Go Ram minimum\n",
"n_geod = 1000\n",
"nx, ny = 720, 360"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"# n_cpu = 36 # 144 Go Ram minimum\n",
"# n_geod = 30000\n",
"# nx, ny = 4000, 2000"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"%display latex"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/latex": [
"\\begin{math}\n",
"\\newcommand{\\Bold}[1]{\\mathbf{#1}}t :\\ \\left( -\\infty, +\\infty \\right) ;\\quad r :\\ \\left( 4.20000000000000 , +\\infty \\right) ;\\quad {\\theta} :\\ \\left( 0 , \\pi \\right) ;\\quad {\\phi} :\\ \\left( -\\infty, +\\infty \\right)\n",
"\\end{math}"
],
"text/plain": [
"t: (-oo, +oo); r: (4.20000000000000, +oo); th: (0, pi); ph: (-oo, +oo)"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"M = Manifold(4, 'M', structure='Lorentzian')\n",
"\n",
"m = var('m')\n",
"a = var('a')\n",
"\n",
"m_val = 2\n",
"a_val = 1\n",
"\n",
"C. = M.chart(r't r:(4.2,+oo) th:(0,pi):\\theta ph:\\phi') #check that minimum radius is bigger than singularity in g[1,1]\n",
"C.coord_range()"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/latex": [
"\\begin{math}\n",
"\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left(\\begin{array}{rrrr}\n",
"\\frac{2 \\, m r}{a^{2} \\cos\\left({\\theta}\\right)^{2} + r^{2}} - 1 & 0 & 0 & -\\frac{2 \\, a m r \\sin\\left({\\theta}\\right)^{2}}{a^{2} \\cos\\left({\\theta}\\right)^{2} + r^{2}} \\\\\n",
"0 & \\frac{a^{2} \\cos\\left({\\theta}\\right)^{2} + r^{2}}{a^{2} - 2 \\, m r + r^{2}} & 0 & 0 \\\\\n",
"0 & 0 & a^{2} \\cos\\left({\\theta}\\right)^{2} + r^{2} & 0 \\\\\n",
"-\\frac{2 \\, a m r \\sin\\left({\\theta}\\right)^{2}}{a^{2} \\cos\\left({\\theta}\\right)^{2} + r^{2}} & 0 & 0 & {\\left(\\frac{2 \\, a^{2} m r \\sin\\left({\\theta}\\right)^{2}}{a^{2} \\cos\\left({\\theta}\\right)^{2} + r^{2}} + a^{2} + r^{2}\\right)} \\sin\\left({\\theta}\\right)^{2}\n",
"\\end{array}\\right)\n",
"\\end{math}"
],
"text/plain": [
"[ 2*m*r/(a^2*cos(th)^2 + r^2) - 1 0 0 -2*a*m*r*sin(th)^2/(a^2*cos(th)^2 + r^2)]\n",
"[ 0 (a^2*cos(th)^2 + r^2)/(a^2 - 2*m*r + r^2) 0 0]\n",
"[ 0 0 a^2*cos(th)^2 + r^2 0]\n",
"[ -2*a*m*r*sin(th)^2/(a^2*cos(th)^2 + r^2) 0 0 (2*a^2*m*r*sin(th)^2/(a^2*cos(th)^2 + r^2) + a^2 + r^2)*sin(th)^2]"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"g = M.metric()\n",
"g[0,0] = (2*m*r)/(r^2+(a*cos(th))^2)-1\n",
"g[0,3] = -2*m*r*a*sin(th)^2/(r^2+(a*cos(th))^2)\n",
"g[1,1] = (r^2+(a*cos(th))^2)/(r^2-2*m*r+a^2)\n",
"g[2,2] = r^2+(a*cos(th))^2\n",
"g[3,3] = (r^2+a^2+(2*m*r*a^2*sin(th)^2)/(r^2+(a*cos(th))^2))*sin(th)^2\n",
"g[:]"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/latex": [
"\\begin{math}\n",
"\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\begin{array}{llcl} & M & \\longrightarrow & \\mathbb{E}^{3} \\\\ & \\left(t, r, {\\theta}, {\\phi}\\right) & \\longmapsto & \\left(x, y, z\\right) = \\left(r \\cos\\left({\\phi}\\right) \\sin\\left({\\theta}\\right), r \\sin\\left({\\phi}\\right) \\sin\\left({\\theta}\\right), r \\cos\\left({\\theta}\\right)\\right) \\end{array}\n",
"\\end{math}"
],
"text/plain": [
"M --> E^3\n",
" (t, r, th, ph) |--> (x, y, z) = (r*cos(ph)*sin(th), r*sin(ph)*sin(th), r*cos(th))"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"E. = EuclideanSpace()\n",
"phi = M.diff_map(E, [r*sin(th)*cos(ph), r*sin(th)*sin(ph), r*cos(th)])\n",
"phi.display()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## One Geodesic"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Try the following for circular Orbit"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n"
],
"text/plain": [
"Graphics3d Object"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"p = M((0, 14.98, pi/2, 0))\n",
"Tp = M.tangent_space(p)\n",
"v = Tp((2, 0, 0.005, 0.05))\n",
"v = v / sqrt(-g.at(p)(v, v))\n",
"\n",
"tau = var('tau')\n",
"\n",
"curve = M.integrated_geodesic(g, (tau, 0, 5000), v)\n",
"sol = curve.solve(step = 1, method=\"ode_int\", parameters_values={m: m_val, a: a_val})\n",
"\n",
"interp = curve.interpolate()\n",
"\n",
"P = curve.plot_integrated(mapping=phi, color=\"red\", thickness=2, plot_points=5000)\n",
"P += sage.plot.plot3d.shapes.Sphere(4, color='grey')\n",
"P"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Try the following for light at different Values for a, aimed directly at the Black Hole"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [],
"source": [
"# p = M((0, 20, pi/2, 0))\n",
"# Tp = M.tangent_space(p)\n",
"# v = Tp((2, -1, 0.000, 0.00))\n",
"# v = v / sqrt(-g.at(p)(v, v))\n",
"\n",
"# A=[0,1,2,4] #list of values for a\n",
"# Color=[[\"red\"], [\"blue\"], [\"green\"],[\"orange\"], [\"black\"]] #color for different Geodesics\n",
"\n",
"# P = sage.plot.plot3d.shapes.Sphere(4, color='grey')\n",
"\n",
"# tau = var('tau')\n",
"# for i in range(4):\n",
"# curve = M.integrated_geodesic(g, (tau, 0, 200), initial_tangent_vector=v, across_charts=True)\n",
"# curve.solve_across_charts(step=0.2, parameters_values={m:m_val, a:A[i]})\n",
"# interp = curve.interpolate()\n",
"# P += curve.plot_integrated(mapping=phi, color=Color[i], thickness=2, plot_points=10000, across_charts=True)\n",
"\n",
"# P"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Multiple Geodesics"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [],
"source": [
"import multiprocessing\n",
"from ipywidgets import FloatProgress\n",
"from IPython.display import display"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [],
"source": [
"def chunks(l, n):\n",
" \"\"\"Yield successive n-sized chunks from l.\"\"\"\n",
" for i in range(0, len(l), n):\n",
" yield l[i:i + n]\n",
"\n",
"n_batches_per_cpu = 3"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
"curve = M.integrated_geodesic(g, (tau, 0, 200), v, across_charts=True)"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [],
"source": [
"args = []\n",
"start_index = 0\n",
"\n",
"for chunk in chunks(range(n_geod), n_geod//(n_batches_per_cpu*n_cpu)):\n",
" args += [(loads(curve.dumps()), start_index, len(chunk))] #\n",
" start_index += len(chunk)"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/latex": [
"\\begin{math}\n",
"\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left[\\mathit{dt} = \\frac{2 \\, {\\left(a^{3} m - 2 \\, a m^{2} r_{0} + a m r_{0}^{2}\\right)} y - \\sqrt{-2 \\, a^{2} m r_{0} - 4 \\, m r_{0}^{3} + r_{0}^{4} + {\\left(a^{2} + 4 \\, m^{2}\\right)} r_{0}^{2} + {\\left(a^{6} - 6 \\, a^{4} m r_{0} - 6 \\, m r_{0}^{5} + r_{0}^{6} + 3 \\, {\\left(a^{2} + 4 \\, m^{2}\\right)} r_{0}^{4} - 4 \\, {\\left(3 \\, a^{2} m + 2 \\, m^{3}\\right)} r_{0}^{3} + 3 \\, {\\left(a^{4} + 4 \\, a^{2} m^{2}\\right)} r_{0}^{2}\\right)} y^{2}} r_{0}}{2 \\, a^{2} m + 4 \\, m r_{0}^{2} - r_{0}^{3} - {\\left(a^{2} + 4 \\, m^{2}\\right)} r_{0}}, \\mathit{dt} = \\frac{2 \\, {\\left(a^{3} m - 2 \\, a m^{2} r_{0} + a m r_{0}^{2}\\right)} y + \\sqrt{-2 \\, a^{2} m r_{0} - 4 \\, m r_{0}^{3} + r_{0}^{4} + {\\left(a^{2} + 4 \\, m^{2}\\right)} r_{0}^{2} + {\\left(a^{6} - 6 \\, a^{4} m r_{0} - 6 \\, m r_{0}^{5} + r_{0}^{6} + 3 \\, {\\left(a^{2} + 4 \\, m^{2}\\right)} r_{0}^{4} - 4 \\, {\\left(3 \\, a^{2} m + 2 \\, m^{3}\\right)} r_{0}^{3} + 3 \\, {\\left(a^{4} + 4 \\, a^{2} m^{2}\\right)} r_{0}^{2}\\right)} y^{2}} r_{0}}{2 \\, a^{2} m + 4 \\, m r_{0}^{2} - r_{0}^{3} - {\\left(a^{2} + 4 \\, m^{2}\\right)} r_{0}}\\right]\n",
"\\end{math}"
],
"text/plain": [
"[dt == (2*(a^3*m - 2*a*m^2*r0 + a*m*r0^2)*y - sqrt(-2*a^2*m*r0 - 4*m*r0^3 + r0^4 + (a^2 + 4*m^2)*r0^2 + (a^6 - 6*a^4*m*r0 - 6*m*r0^5 + r0^6 + 3*(a^2 + 4*m^2)*r0^4 - 4*(3*a^2*m + 2*m^3)*r0^3 + 3*(a^4 + 4*a^2*m^2)*r0^2)*y^2)*r0)/(2*a^2*m + 4*m*r0^2 - r0^3 - (a^2 + 4*m^2)*r0), dt == (2*(a^3*m - 2*a*m^2*r0 + a*m*r0^2)*y + sqrt(-2*a^2*m*r0 - 4*m*r0^3 + r0^4 + (a^2 + 4*m^2)*r0^2 + (a^6 - 6*a^4*m*r0 - 6*m*r0^5 + r0^6 + 3*(a^2 + 4*m^2)*r0^4 - 4*(3*a^2*m + 2*m^3)*r0^3 + 3*(a^4 + 4*a^2*m^2)*r0^2)*y^2)*r0)/(2*a^2*m + 4*m*r0^2 - r0^3 - (a^2 + 4*m^2)*r0)]"
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"dt, y, r0 = var('dt, y, r0')\n",
"\n",
"p = M((0, r0, pi/2, 0))\n",
"Tp = M.tangent_space(p)\n",
"v = Tp((dt, -1, 0, y))\n",
"\n",
"sol = g.at(p)(v, v).solve(dt)\n",
"sol"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [],
"source": [
"def calc_some_geodesics(args):\n",
" \"\"\"\n",
" Compute nb geodesics starting at index n0\n",
" \"\"\"\n",
" curve, n0, nb = args\n",
" res = {}\n",
" r = 100\n",
" posi = [0, r, pi/2, 0]\n",
" p = M(posi)\n",
" Tp = M.tangent_space(p)\n",
" for i in range(n0, n0+nb):\n",
" dy = i*0.006/n_geod\n",
" v = Tp([sol[0].rhs()(r0=r, y=dy, m=m_val, a=a_val).n(), -1, 0, dy])\n",
" curve._initial_tangent_vector = v\n",
" curve.solve_across_charts(step=0.2, parameters_values={m:m_val, a:a_val})\n",
" res[i] = (p.coord(), curve._solutions.copy())\n",
" curve._solutions.clear()\n",
" return res"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "fef6a25ffd99454da6f5a0bada40f10b",
"version_major": 2,
"version_minor": 0
},
"text/plain": [
"FloatProgress(value=0.0, max=1000.0)"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"geo = {}\n",
"pool = multiprocessing.Pool(n_cpu)\n",
"\n",
"%display plain\n",
"f = FloatProgress(min=0, max=n_geod)\n",
"display(f)\n",
"\n",
"for i, some_res in enumerate(pool.map(calc_some_geodesics, args)):\n",
" f.value += len(some_res)\n",
" geo.update(some_res)\n",
"\n",
"pool.close()\n",
"pool.join()"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n"
],
"text/plain": [
"Graphics3d Object"
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"P = sage.plot.plot3d.shapes.Sphere(4, color='grey')\n",
"\n",
"for i in range(0, n_geod, 5*n_geod/100): \n",
" curve._solutions = geo[i][1]\n",
" interp = curve.interpolate()\n",
" P += curve.plot_integrated(mapping=phi, color=[\"red\"], thickness=2, plot_points=150, \n",
" label_axes=False, across_charts=True)#\n",
"\n",
"P"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"8.46600505906165\n"
]
}
],
"source": [
"Angle = pi/20 #angle used in the actual calculation\n",
"\n",
"Xi=a_val/m_val\n",
"Z_1=1+(1-Xi^2)^(1/3)*((1+Xi)^(1/3)+(1-Xi)^(1/3))\n",
"Z_2=(3*Xi^2+Z_1^2)^(1/2)\n",
"#calculating ISCO\n",
"\n",
"disk_min = m_val*(3+Z_2-((3-Z_1)*(3+Z_1+2*Z_2))^(1/2))\n",
"print(numerical_approx(disk_min))\n",
"if disk_min<2*m_val:\n",
" disk_min=2*m_val\n",
"\n",
"disk_max = (50^2+disk_min^2-12^2)^(1/2)\n",
"#disk_max =25\n",
"#keeping area of accretion Disk the same as in original work of Florentin Jaffredo\n",
"\n",
"#disk_min = 4\n",
"#disk_max = 50\n",
"alpha = -pi/20\n"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {
"scrolled": true
},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n"
],
"text/plain": [
"Graphics3d Object"
]
},
"execution_count": 19,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"D = sage.plot.plot3d.shapes.Torus((disk_min+disk_max)/2,\n",
" (disk_max-disk_min)/2).scale(1,1,0.01).rotateY(-pi/20)\n",
"\n",
"P + D"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n"
],
"text/plain": [
"Graphics3d Object"
]
},
"execution_count": 20,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"P + D.rotateX(pi/3)"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [],
"source": [
"#now begins the part where we calculate intersection with accretion disk and visualize the image\n",
"geo = [list(geo[i][1].values())[0][0][1].tolist() for i in range(len(geo))]"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [],
"source": [
"def intersection(curve, alpha, beta):\n",
" \"\"\"\n",
" Return True if the curve intersect the disk comprised between dmin and dmax\n",
" tilted of angles alpha and beta\n",
" \"\"\"\n",
" n = len(curve)\n",
" r, theta, phi = curve[0][2:5]\n",
" x, y, z = r*sin(theta)*cos(phi), r*sin(theta)*sin(phi), r*cos(theta)\n",
" x, y, z = x, y*cos(beta)-z*sin(beta), z*cos(beta)+y*sin(beta)\n",
" z = z*cos(alpha)+x*sin(alpha)\n",
" for i in range(1, n):\n",
" #Â done in 3 lines for speed consideration\n",
" r = curve[i][2]\n",
" theta = curve[i][3]\n",
" phi = curve[i][4]\n",
" # conversion to cartesian:\n",
" x2, y2, z2 = r*sin(theta)*cos(phi), r*sin(theta)*sin(phi), r*cos(theta) \n",
" # rotation around the X-axis:\n",
" y2, z2 = y2*cos(beta)-z2*sin(beta), z2*cos(beta)+y2*sin(beta) \n",
" # rotation around the Y-axis:\n",
" x2, z2 = x2*cos(alpha)-z2*sin(alpha), z2*cos(alpha)+x2*sin(alpha) \n",
" if z!=z2: # needed to prevent a division by zero next line\n",
" t = z/(z-z2) # if 0<=t<1 then the curve intersect the disk between the points i and i-1\n",
" if t>=0 and t<1 and curve[i][2]>disk_min and curve[i][2]=0 and t<1 and curve[i][2]>dmin and curve[i][2]>> from scipy.misc import bytescale\n",
" >>> img = np.array([[ 91.06794177, 3.39058326, 84.4221549 ],\n",
" ... [ 73.88003259, 80.91433048, 4.88878881],\n",
" ... [ 51.53875334, 34.45808177, 27.5873488 ]])\n",
" >>> bytescale(img)\n",
" array([[255, 0, 236],\n",
" [205, 225, 4],\n",
" [140, 90, 70]], dtype=uint8)\n",
" >>> bytescale(img, high=200, low=100)\n",
" array([[200, 100, 192],\n",
" [180, 188, 102],\n",
" [155, 135, 128]], dtype=uint8)\n",
" >>> bytescale(img, cmin=0, cmax=255)\n",
" array([[91, 3, 84],\n",
" [74, 81, 5],\n",
" [52, 34, 28]], dtype=uint8)\n",
" \"\"\"\n",
" if data.dtype == np.uint8:\n",
" return data\n",
"\n",
" if high > 255:\n",
" raise ValueError(\"`high` should be less than or equal to 255.\")\n",
" if low < 0:\n",
" raise ValueError(\"`low` should be greater than or equal to 0.\")\n",
" if high < low:\n",
" raise ValueError(\"`high` should be greater than or equal to `low`.\")\n",
"\n",
" if cmin is None:\n",
" cmin = data.min()\n",
" if cmax is None:\n",
" cmax = data.max()\n",
"\n",
" cscale = cmax - cmin\n",
" if cscale < 0:\n",
" raise ValueError(\"`cmax` should be larger than `cmin`.\")\n",
" elif cscale == 0:\n",
" cscale = 1\n",
"\n",
" scale = float(high - low) / cscale\n",
" bytedata = (data - cmin) * scale + low\n",
" return (bytedata.clip(low, high) + 0.5).astype(np.uint8)\n",
"\n",
"def toimage(arr, high=255, low=0, cmin=None, cmax=None, pal=None,\n",
" mode=None, channel_axis=None):\n",
" \"\"\"Takes a numpy array and returns a PIL image.\n",
" This function is only available if Python Imaging Library (PIL) is installed.\n",
" The mode of the PIL image depends on the array shape and the `pal` and\n",
" `mode` keywords.\n",
" For 2-D arrays, if `pal` is a valid (N,3) byte-array giving the RGB values\n",
" (from 0 to 255) then ``mode='P'``, otherwise ``mode='L'``, unless mode\n",
" is given as 'F' or 'I' in which case a float and/or integer array is made.\n",
" .. warning::\n",
" This function uses `bytescale` under the hood to rescale images to use\n",
" the full (0, 255) range if ``mode`` is one of ``None, 'L', 'P', 'l'``.\n",
" It will also cast data for 2-D images to ``uint32`` for ``mode=None``\n",
" (which is the default).\n",
" Notes\n",
" -----\n",
" For 3-D arrays, the `channel_axis` argument tells which dimension of the\n",
" array holds the channel data.\n",
" For 3-D arrays if one of the dimensions is 3, the mode is 'RGB'\n",
" by default or 'YCbCr' if selected.\n",
" The numpy array must be either 2 dimensional or 3 dimensional.\n",
" \"\"\"\n",
" data = np.asarray(arr)\n",
" if np.iscomplexobj(data):\n",
" raise ValueError(\"Cannot convert a complex-valued array.\")\n",
" shape = list(data.shape)\n",
" valid = len(shape) == 2 or ((len(shape) == 3) and\n",
" ((3 in shape) or (4 in shape)))\n",
" if not valid:\n",
" raise ValueError(\"'arr' does not have a suitable array shape for \"\n",
" \"any mode.\")\n",
" if len(shape) == 2:\n",
" shape = (shape[1], shape[0]) # columns show up first\n",
" if mode == 'F':\n",
" data32 = data.astype(np.float32)\n",
" image = Image.frombytes(mode, shape, data32.tobytes())\n",
" return image\n",
" if mode in [None, 'L', 'P']:\n",
" bytedata = bytescale(data, high=high, low=low,\n",
" cmin=cmin, cmax=cmax)\n",
" image = Image.frombytes('L', shape, bytedata.tobytes())\n",
" if pal is not None:\n",
" image.putpalette(np.asarray(pal, dtype=np.uint8).tobytes())\n",
" # Becomes a mode='P' automagically.\n",
" elif mode == 'P': # default gray-scale\n",
" pal = (np.arange(0, 256, 1, dtype=np.uint8)[:, np.newaxis] *\n",
" np.ones((3,), dtype=np.uint8)[np.newaxis, :])\n",
" image.putpalette(np.asarray(pal, dtype=np.uint8).tobytes())\n",
" return image\n",
" if mode == '1': # high input gives threshold for 1\n",
" bytedata = (data > high)\n",
" image = Image.frombytes('1', shape, bytedata.tobytes())\n",
" return image\n",
" if cmin is None:\n",
" cmin = np.amin(np.ravel(data))\n",
" if cmax is None:\n",
" cmax = np.amax(np.ravel(data))\n",
" data = (data*1.0 - cmin)*(high - low)/(cmax - cmin) + low\n",
" if mode == 'I':\n",
" data32 = data.astype(np.uint32)\n",
" image = Image.frombytes(mode, shape, data32.tobytes())\n",
" else:\n",
" raise ValueError(_errstr)\n",
" return image\n",
"\n",
" # if here then 3-d array with a 3 or a 4 in the shape length.\n",
" # Check for 3 in datacube shape --- 'RGB' or 'YCbCr'\n",
" if channel_axis is None:\n",
" if (3 in shape):\n",
" ca = np.flatnonzero(np.asarray(shape) == 3)[0]\n",
" else:\n",
" ca = np.flatnonzero(np.asarray(shape) == 4)\n",
" if len(ca):\n",
" ca = ca[0]\n",
" else:\n",
" raise ValueError(\"Could not find channel dimension.\")\n",
" else:\n",
" ca = channel_axis\n",
"\n",
" numch = shape[ca]\n",
" if numch not in [3, 4]:\n",
" raise ValueError(\"Channel axis dimension is not valid.\")\n",
"\n",
" bytedata = bytescale(data, high=high, low=low, cmin=cmin, cmax=cmax)\n",
" if ca == 2:\n",
" strdata = bytedata.tobytes()\n",
" shape = (shape[1], shape[0])\n",
" elif ca == 1:\n",
" strdata = np.transpose(bytedata, (0, 2, 1)).tobytes()\n",
" shape = (shape[2], shape[0])\n",
" elif ca == 0:\n",
" strdata = np.transpose(bytedata, (1, 2, 0)).tobytes()\n",
" shape = (shape[2], shape[1])\n",
" if mode is None:\n",
" if numch == 3:\n",
" mode = 'RGB'\n",
" else:\n",
" mode = 'RGBA'\n",
"\n",
" if mode not in ['RGB', 'RGBA', 'YCbCr', 'CMYK']:\n",
" raise ValueError(_errstr)\n",
"\n",
" if mode in ['RGB', 'YCbCr']:\n",
" if numch != 3:\n",
" raise ValueError(\"Invalid array shape for mode.\")\n",
" if mode in ['RGBA', 'CMYK']:\n",
" if numch != 4:\n",
" raise ValueError(\"Invalid array shape for mode.\")\n",
"\n",
" # Here we know data and mode is correct\n",
" image = Image.frombytes(mode, shape, strdata)\n",
" return image"
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAtAAAAFoCAIAAADxRFtOAAAMXklEQVR4nO3d23bbOBIFUCur//+XPQ9OPLKsC0WygKrC3g9ZmVndjkSiDo5Bpf3xAQAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAMBGl9kvAMjo89i/LlmAG2IBFnWwUhwhd2BBBh+am1gs3iWPoDEDDg0VKhmPyCZoxlBDeQ3qxRbSCkozwlDVIj3jN7EFFZlcKGPZhvGcFIMSjCpkp2dsJM4gsz+zXwBwx+e/X7WN7T6vrhuQjW8JIBeb5YkEHORhHiELVSOImIMMTCJMpmcMI+9gIgMIc+gZEwk+GM/cwWiqRhLiD0YycTCInpGWHIQBDBqEUzVKkIYQyohBIFWjHJkIQQwXhFA1SpOMcDpjBSdTNdqQj3AiAwWnUTVakpJwCqMEJ1A12pOVcJAhgkNUjaVITNjN+MBOqsay5Cbs4MfTwx7axsrcfdhBU4f32Gz4JkBhO/MCW6ka3CVGYQuPVGATbYNHrA3YQjWHF2wnbCRP4QkDAg+pGuwgVeEuj1TgPm2DfawcuEsXh1s2DE4hXuGaEw74QdvgLNYSXFPB4S/bA0HkLHwYBPhQNRhC2rI4j1RYnbbBGFYai1M4WJo9gJGsN1bmkI9FiX4mkrwsyAkHK9I2mMsKZEEKB8uR9WRgHbIaB3ssRMSTkBRmEU44WIW2QU5WJotQOFiCTCcz65MVOMyjOVFOIRKZxpxw0Jm2QS1WLI3p0/QkuClNNNOPEw4a0jaozhqmH4WDbiQ1PVjJNOPcjj4ENC2JaXpwwkET2gZdWdv0oHDQgUSmNyucBhQOypPFrMA6pzoPBylMBLMgqU1RTjgAgHC6MiU522BxsptynHBQj7YBpoBytGQqEbJwQ4hThRMOytA24DdzQRUKBzVIVXjEdFCCwkEB8hSeMyPk5/EfqYlReItMJy0nHOSlbcC7TA1pKRwkJTdhH7NDTgoHGUlMOMIEkZDCQTqyEo4zR2TjA0YkIiLhdFKeJJxwkIW2ARFMFkkoHKQgEyGO+SIDhYP5pCFEM2VM5+keMwlBGEzoM4sTDqbRNmA8c8csCgdzSD2YxfQxhcLBBPIO5jKDjOdxHkOJOUjFHsAwFhuDqBqQlp2AATxSAQDC6bWEc7ZRzim54L6XYz8glAVGLLtOCQOCwEoowZZAHKuLKDaY5CYOv7WRnI2BCNYVIewoaaWaeeskrVTrhB4sKk5mC0kr7bRbM2mlXTNUZDlxJjtHQoWG3PpJqND6ITlriRPYJ3IqOt6WU05FlxN5WEIcZXtIqMFgW1cJNVhXTGT9sJ8tIaFmI22NJdRsjTGMlcNOdoJsGg+zxZZN48VGHMuGt0n/hNpPslWXUPtVx7ksGN4g9BOKm+HdtzvhSyKOXYSNLBU2EfQ5nTvAQXe5xIvkIHsJL1kkvCDf0zpreh/d4t1f//QvuPHrM50dhScsD56R7DkF/TTXoDgI+oMszpxsKjxibXCHKM/s4NAO6xkD/mgLNTO7CzcsCX6Q4Mkdmdibm5vnp8We+KbIxh7DN4uBvwR3fmd9qCLJ2J/1qizd/JIsOeayDJDXNeyb1eubm3baj79Ia7iEtCuQMSyAdcnoQhq3jS86x1KSr0aCuO8rEs2FtK8a19SOpRRamZzCHV+LOK7l4KZbdLwPvn6LvJaiq5Qd3OsliOCKdgxn0YON3w6+EQu+otIrli3c4rZkbmlHdtk2U33kHVn/pbVZw1xzWxsStdVpG990jpU1W8y4oa1I2B7eGss2j1Ee2f0GjUMPLVf1mtzK8qRqM/v21PaTvO+dmo5m2q/z3ty+qiRpS9rGEzoH3xZZ8824a8VIz8a0jZd0Dm4stf6rc7NqkJjtaRsb6RzcteAslOMe5SUi16FtvEXn4LllRyM59yUdsbig7XOobXzZcR1M1oIWH5Ns3I755ODitI19dA7eZXDm+m/2C1iX7AMYSWWfy2UfSsnghuONIxxycJyBGsalDifgeETbOE7n4ETmK5TLG0Ki8ZK2cRadgwjG7XQ+w3EaEQbQRvufUjSey3iIksE+jjfO5ZCDYUzibi7dHqKKgzYOnrax3bvXyhRzkKl8l0cqW4knziKnMrgYao7xzOVdrtJDwoggjjeCOOQgCTN7lxOOHwQQ0SRRHg45COLw4y6XQuIwlOONUA45yGzxcV707UsZptA2BtA5KGHB6V7okYpYASCJBR+7NH+bSgZ5vPufiGg+nJHevYCCgjwaD37DEw7ZQWkW8Fk+W2c3XTU++WhSOGQ0zTQLmsH89RN6aFY+Cr8FgUIhPi46mI+O0ljRfCh2wiEUAFhc0ZOPAi9VyaA6xxtTOORgNcmjI+8Jh+EHgO2Sf9OSq3AoGfSTc/L5zUdNaSPnM5f5r8SE05vnKRN5qgLX5sbLtBMOgw0AI839xmZo4VAyWI3jilo8VWERU565jPiDDDDL8jxlOk9VYKPo/Ak84TC3AFBF9Hc+JxcOJQPeZWrG8KNVYKOgBy7nFA6JCb/Z3iryMQ64duKxx6GvYCzhiS3TlfOvyzfz7kWWbPDE7qTac8JhGuEURmk8D1bgoN1nHu8VDvkIAHz8qwTba8emf1LPgB2eT9fNWPm2O9RbV1viwQ4vQ+zZCYepgyOc3lck92Cfl49a/jz610wdHPdojszXXO4LxHlUIS43/xBwrt9l//egOQgZYMtll4Fwusv1b8wYRPseudfFnxjPr7wYhGh/jBkM8Hn1K3m4LzDM/c9wAKfbsqvZ+c7lmkMeCgdM5nnKGK4zzKVwAADhAn88PbDFu/+1Pnbw3ASmu5hDSEj/OEiyQTYKB2SnfGwkzSAzhQOqWraISC2oSOGAtoIaye4fTv3u1wc6UTgAgHD+WiwAEE7hAADCKRwAQDiFAwAIp3AAAOEUDgAgnMIBAIRTOACAcAoHABBO4QAAwikcAEA4hQMACKdwAADhFA4AIJzCAQCEUzgAgHAKBwAQTuEAAMIpHABAOIUDAAincAAA4RQOACCcwgEAhFM4AIBwCgcAEE7hAADCKRwAQDiFAwAIp3AAAOEUDgAgnMIBAIRTOACAcAoHABBO4QAAwikcAEA4hQMACKdwAADhFA4AIJzCAQCEUzgAgHAKBwAQTuEAAMIpHABAOIUDAAincAAA4RQOACCcwgEAhFM4AIBwCgcAEE7hAADCXb5/9znxVQAATV2ufv1B8wAADro8/Z8/aB4AwFseFYtnheOL2gEAvPS8UrwuHN80DwDgxsYm8Ubh+KZ5AMDi3i0QewrHNeUDABZxpDQcLRxf1A4AaOx4XTincFxTPgCggXMrwvmF44b+AQAlhHaC8MLxRe0AgLQGtIFBheOG/gEAE43f/ucUjhv6BwCEmr7fT38B/6d2AMDpkuz0SV7GfSoIALwl7b6e9oX9oHkAwBP5t/P8r/A+FQSAZVXcvCu+5h8+Pz4u+gcArV3+7Xd1lX7x9ykfADTQbIdu9nZuOf8AoIQGZxjPNX5rD+kfAEy32ga82vu9pXwAMMzKm+7K7/0ZRQSA3Wyuv7kmm+gfADxhN33JJTpKFwFYhC3zCFfvTMoHQDO2ybO4kkNpJACp2AWHcaln0j8ABrPtzeLK56KCAJzIJpeHe1GPUgLwYQOrxv3qQAUB2rNdVecOtqWFAEXZmVpyW/lLQQGC2Gn4sAx4SREBNrKj8ITlwX66CCzItsE+Vg6BvhrJ5d7/CUx3dzbtCgSxtEhKL4HdJDsJWZYUppSwIKlNUZYuq9BOSEsQswLrHF74/DcnN7/5fuD9ee9XqntyZy+PVwXwiAGBjFSWg0QbZPNn9gsA7rBfHuHqQUIGE1Jz1PEWiQZpGU8oQO14SZZBcoYUylA77pJiUIJRhWLUjm/yCwrxoVEoxi77xXWAWswsFLbgaYfMgqIML3TQvnmIKqjOIxXooPd+3PvdwSIMMjTU4MBDNkEzhho6K9c8RBJ0ZbphFWnLhxiCFZh0WE6S5iF9YClGHhhRQWQNLE4IAC9srCPSBAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAADI5H+NqLVGCgyWmAAAAABJRU5ErkJggg==\n",
"text/plain": [
""
]
},
"execution_count": 31,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"img1 = toimage(data)\n",
"img1"
]
},
{
"cell_type": "code",
"execution_count": 32,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/latex": [
"\\begin{math}\n",
"\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left(\\begin{array}{rrrr}\n",
"\\frac{2 \\, m r}{a^{2} \\cos\\left({\\theta}\\right)^{2} + r^{2}} - 1 & 0 & 0 & -\\frac{2 \\, a m r \\sin\\left({\\theta}\\right)^{2}}{a^{2} \\cos\\left({\\theta}\\right)^{2} + r^{2}} \\\\\n",
"0 & \\frac{a^{2} \\cos\\left({\\theta}\\right)^{2} + r^{2}}{a^{2} - 2 \\, m r + r^{2}} & 0 & 0 \\\\\n",
"0 & 0 & a^{2} \\cos\\left({\\theta}\\right)^{2} + r^{2} & 0 \\\\\n",
"-\\frac{2 \\, a m r \\sin\\left({\\theta}\\right)^{2}}{a^{2} \\cos\\left({\\theta}\\right)^{2} + r^{2}} & 0 & 0 & {\\left(\\frac{2 \\, a^{2} m r \\sin\\left({\\theta}\\right)^{2}}{a^{2} \\cos\\left({\\theta}\\right)^{2} + r^{2}} + a^{2} + r^{2}\\right)} \\sin\\left({\\theta}\\right)^{2}\n",
"\\end{array}\\right)\n",
"\\end{math}"
],
"text/plain": [
"[ 2*m*r/(a^2*cos(th)^2 + r^2) - 1 0 0 -2*a*m*r*sin(th)^2/(a^2*cos(th)^2 + r^2)]\n",
"[ 0 (a^2*cos(th)^2 + r^2)/(a^2 - 2*m*r + r^2) 0 0]\n",
"[ 0 0 a^2*cos(th)^2 + r^2 0]\n",
"[ -2*a*m*r*sin(th)^2/(a^2*cos(th)^2 + r^2) 0 0 (2*a^2*m*r*sin(th)^2/(a^2*cos(th)^2 + r^2) + a^2 + r^2)*sin(th)^2]"
]
},
"execution_count": 32,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"%display latex\n",
"g[:]"
]
},
{
"cell_type": "code",
"execution_count": 33,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/latex": [
"\\begin{math}\n",
"\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left(\\begin{array}{rrrr}\n",
"\\frac{2 \\, m - r_{0}}{r_{0}} & 0 & 0 & -\\frac{2 \\, a m}{r_{0}} \\\\\n",
"0 & \\frac{r_{0}^{2}}{a^{2} - 2 \\, m r_{0} + r_{0}^{2}} & 0 & 0 \\\\\n",
"0 & 0 & r_{0}^{2} & 0 \\\\\n",
"-\\frac{2 \\, a m}{r_{0}} & 0 & 0 & \\frac{2 \\, a^{2} m + a^{2} r_{0} + r_{0}^{3}}{r_{0}}\n",
"\\end{array}\\right)\n",
"\\end{math}"
],
"text/plain": [
"[ (2*m - r_0)/r_0 0 0 -2*a*m/r_0]\n",
"[ 0 r_0^2/(a^2 - 2*m*r_0 + r_0^2) 0 0]\n",
"[ 0 0 r_0^2 0]\n",
"[ -2*a*m/r_0 0 0 (2*a^2*m + a^2*r_0 + r_0^3)/r_0]"
]
},
"execution_count": 33,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"r0, phi = var('r_0, phi')\n",
"p = M((0, r0, pi/2, phi))\n",
"g.at(p)[:]"
]
},
{
"cell_type": "code",
"execution_count": 34,
"metadata": {},
"outputs": [],
"source": [
"# default frame\n",
"fr = C.frame()\n",
"\n",
"# create an automorphism field\n",
"aut = M.automorphism_field()\n",
"\n",
"# some symbolic variables\n",
"a, b, c = var('a, b, c')\n",
"\n",
"# let's try with the simplest matrix possible\n",
"aut.add_comp()[:] = [[a, 0, 0, 0], [b, c, 0, 0], [0, 0, 1/r0, 0], \n",
" [0, 0, 0, 1/r0]] # only b is off-diagonal\n",
"fr2 = fr.new_frame(aut, 'f2')"
]
},
{
"cell_type": "code",
"execution_count": 35,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/latex": [
"\\begin{math}\n",
"\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left(\\begin{array}{rrrr}\n",
"\\frac{b^{2} r_{0}^{2}}{a^{2} - 2 \\, m r_{0} + r_{0}^{2}} + \\frac{a^{2} {\\left(2 \\, m - r_{0}\\right)}}{r_{0}} & \\frac{b c r_{0}^{2}}{a^{2} - 2 \\, m r_{0} + r_{0}^{2}} & 0 & -\\frac{2 \\, a^{2} m}{r_{0}^{2}} \\\\\n",
"\\frac{b c r_{0}^{2}}{a^{2} - 2 \\, m r_{0} + r_{0}^{2}} & \\frac{c^{2} r_{0}^{2}}{a^{2} - 2 \\, m r_{0} + r_{0}^{2}} & 0 & 0 \\\\\n",
"0 & 0 & 1 & 0 \\\\\n",
"-\\frac{2 \\, a^{2} m}{r_{0}^{2}} & 0 & 0 & \\frac{2 \\, a^{2} m + a^{2} r_{0} + r_{0}^{3}}{r_{0}^{3}}\n",
"\\end{array}\\right)\n",
"\\end{math}"
],
"text/plain": [
"[b^2*r_0^2/(a^2 - 2*m*r_0 + r_0^2) + a^2*(2*m - r_0)/r_0 b*c*r_0^2/(a^2 - 2*m*r_0 + r_0^2) 0 -2*a^2*m/r_0^2]\n",
"[ b*c*r_0^2/(a^2 - 2*m*r_0 + r_0^2) c^2*r_0^2/(a^2 - 2*m*r_0 + r_0^2) 0 0]\n",
"[ 0 0 1 0]\n",
"[ -2*a^2*m/r_0^2 0 0 (2*a^2*m + a^2*r_0 + r_0^3)/r_0^3]"
]
},
"execution_count": 35,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"g.at(p)[fr2.at(p), :]"
]
},
{
"cell_type": "code",
"execution_count": 36,
"metadata": {},
"outputs": [],
"source": [
"c = sqrt(r/(r+2*m))\n",
"a = sqrt(((r+2*m)/(r)))\n",
"b = -2*a*m/(2*m+r)"
]
},
{
"cell_type": "code",
"execution_count": 37,
"metadata": {},
"outputs": [],
"source": [
"aut2 = M.automorphism_field() # new automorphism field\n",
"aut2.add_comp()[:] = [[a, 0, 0, 0], [b, c, 0, 0], [0, 0, 1/r, 0], [0, 0, 0, 1/(r*sin(th))]]\n",
"fr3 = fr.new_frame(aut2, 'f3')"
]
},
{
"cell_type": "code",
"execution_count": 38,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/latex": [
"\\begin{math}\n",
"\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left(\\begin{array}{rrrr}\n",
"\\frac{4 \\, a^{4} m^{2} \\cos\\left({\\theta}\\right)^{4} + 8 \\, a^{2} m^{3} r - 2 \\, a^{2} m r^{3} - r^{6} - {\\left(a^{2} - 12 \\, m^{2}\\right)} r^{4} + 4 \\, {\\left(a^{2} m^{2} - 4 \\, m^{4}\\right)} r^{2} - {\\left(4 \\, a^{4} m^{2} + 2 \\, a^{2} m r^{3} + a^{2} r^{4} + {\\left(a^{4} - 12 \\, a^{2} m^{2}\\right)} r^{2} + 4 \\, {\\left(a^{4} m - 2 \\, a^{2} m^{3}\\right)} r\\right)} \\cos\\left({\\theta}\\right)^{2}}{2 \\, a^{2} m r^{3} + r^{6} + {\\left(a^{2} - 4 \\, m^{2}\\right)} r^{4} + {\\left(2 \\, a^{4} m r + a^{2} r^{4} + {\\left(a^{4} - 4 \\, a^{2} m^{2}\\right)} r^{2}\\right)} \\cos\\left({\\theta}\\right)^{2}} & -\\frac{2 \\, {\\left(a^{2} m \\cos\\left({\\theta}\\right)^{2} + m r^{2}\\right)}}{2 \\, a^{2} m + r^{3} + {\\left(a^{2} - 4 \\, m^{2}\\right)} r} & 0 & -\\frac{2 \\, a \\sqrt{2 \\, m + r} m \\sin\\left({\\theta}\\right)}{{\\left(a^{2} \\cos\\left({\\theta}\\right)^{2} + r^{2}\\right)} \\sqrt{r}} \\\\\n",
"-\\frac{2 \\, {\\left(a^{2} m \\cos\\left({\\theta}\\right)^{2} + m r^{2}\\right)}}{2 \\, a^{2} m + r^{3} + {\\left(a^{2} - 4 \\, m^{2}\\right)} r} & \\frac{a^{2} r \\cos\\left({\\theta}\\right)^{2} + r^{3}}{2 \\, a^{2} m + r^{3} + {\\left(a^{2} - 4 \\, m^{2}\\right)} r} & 0 & 0 \\\\\n",
"0 & 0 & \\frac{a^{2} \\cos\\left({\\theta}\\right)^{2} + r^{2}}{r^{2}} & 0 \\\\\n",
"-\\frac{2 \\, a \\sqrt{2 \\, m + r} m \\sin\\left({\\theta}\\right)}{{\\left(a^{2} \\cos\\left({\\theta}\\right)^{2} + r^{2}\\right)} \\sqrt{r}} & 0 & 0 & \\frac{2 \\, a^{2} m r \\sin\\left({\\theta}\\right)^{2} + a^{2} r^{2} + r^{4} + {\\left(a^{4} + a^{2} r^{2}\\right)} \\cos\\left({\\theta}\\right)^{2}}{a^{2} r^{2} \\cos\\left({\\theta}\\right)^{2} + r^{4}}\n",
"\\end{array}\\right)\n",
"\\end{math}"
],
"text/plain": [
"[(4*a^4*m^2*cos(th)^4 + 8*a^2*m^3*r - 2*a^2*m*r^3 - r^6 - (a^2 - 12*m^2)*r^4 + 4*(a^2*m^2 - 4*m^4)*r^2 - (4*a^4*m^2 + 2*a^2*m*r^3 + a^2*r^4 + (a^4 - 12*a^2*m^2)*r^2 + 4*(a^4*m - 2*a^2*m^3)*r)*cos(th)^2)/(2*a^2*m*r^3 + r^6 + (a^2 - 4*m^2)*r^4 + (2*a^4*m*r + a^2*r^4 + (a^4 - 4*a^2*m^2)*r^2)*cos(th)^2) -2*(a^2*m*cos(th)^2 + m*r^2)/(2*a^2*m + r^3 + (a^2 - 4*m^2)*r) 0 -2*a*sqrt(2*m + r)*m*sin(th)/((a^2*cos(th)^2 + r^2)*sqrt(r))]\n",
"[ -2*(a^2*m*cos(th)^2 + m*r^2)/(2*a^2*m + r^3 + (a^2 - 4*m^2)*r) (a^2*r*cos(th)^2 + r^3)/(2*a^2*m + r^3 + (a^2 - 4*m^2)*r) 0 0]\n",
"[ 0 0 (a^2*cos(th)^2 + r^2)/r^2 0]\n",
"[ -2*a*sqrt(2*m + r)*m*sin(th)/((a^2*cos(th)^2 + r^2)*sqrt(r)) 0 0 (2*a^2*m*r*sin(th)^2 + a^2*r^2 + r^4 + (a^4 + a^2*r^2)*cos(th)^2)/(a^2*r^2*cos(th)^2 + r^4)]"
]
},
"execution_count": 38,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"g[fr3, :]"
]
},
{
"cell_type": "code",
"execution_count": 39,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/latex": [
"\\begin{math}\n",
"\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left(\\begin{array}{rrrr}\n",
"\\frac{\\sqrt{r}}{\\sqrt{2 \\, m + r}} & 0 & 0 & 0 \\\\\n",
"\\frac{2 \\, m}{\\sqrt{2 \\, m + r} \\sqrt{r}} & \\frac{\\sqrt{2 \\, m + r}}{\\sqrt{r}} & 0 & 0 \\\\\n",
"0 & 0 & r & 0 \\\\\n",
"0 & 0 & 0 & r \\sin\\left({\\theta}\\right)\n",
"\\end{array}\\right)\n",
"\\end{math}"
],
"text/plain": [
"[ sqrt(r)/sqrt(2*m + r) 0 0 0]\n",
"[2*m/(sqrt(2*m + r)*sqrt(r)) sqrt(2*m + r)/sqrt(r) 0 0]\n",
"[ 0 0 r 0]\n",
"[ 0 0 0 r*sin(th)]"
]
},
"execution_count": 39,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"aut2.inverse()[:]"
]
},
{
"cell_type": "code",
"execution_count": 40,
"metadata": {},
"outputs": [],
"source": [
"%%cython\n",
"from libc.math cimport cos, sin, sqrt\n",
"\n",
"cpdef tuple spherical_to_xyz(float dr, float dtheta, float dphi, float r, \n",
" float theta, float phi, float alpha, float beta):\n",
" \"\"\"\n",
" Convert spherical coordinates to cartesian and apply the \n",
" two rotations at the same time.\n",
" \"\"\"\n",
" cdef float dx, dy, dz\n",
" cdef float ca, cb, ct, cp\n",
" cdef float sa, sb, st, sp\n",
" \n",
" ca = cos(alpha); sa = sin(alpha)\n",
" cb = cos(beta); sb = sin(beta)\n",
" ct = cos(theta); st = sin(theta)\n",
" cp = cos(phi); sp = sin(phi)\n",
" \n",
" dx = ((-cb*ct*sa - (sa*sb*sp - ca*cp)*st)*dr + \n",
" (r*cb*sa*st - (sa*sb*sp - ca*cp)*r*ct)*dtheta +\n",
" (-(cp*sa*sb + ca*sp)*r*st)*dphi)\n",
" \n",
" dy = ((cb*sp*st - sb*ct)*dr +\n",
" (r*ct*cb*sp + r*sb*st)*dtheta +\n",
" (r*cp*cb*st)*dphi)\n",
" \n",
" dz = ((ca*cb*ct + (ca*sb*sp + cp*sa)*st)*dr +\n",
" (-r*ca*cb*st+(ca*sb*sp + cp*sa)*r*ct)*dtheta +\n",
" ((ca*cp*sb - sa*sp)*r*st)*dphi)\n",
" \n",
" return (dx, dy, dz)\n",
"\n",
"\n",
"cpdef tuple xyz_to_spherical(float dx, float dy, float dz, float x, \n",
" float y, float z):\n",
" \"\"\"\n",
" Convert cartesian back to spherical\n",
" \"\"\"\n",
" cdef r, dr, dth, dph\n",
" r = sqrt(x**2+y**2+z**2)\n",
" dr = (x*dx+y*dy*z*dz)/r\n",
" dth = ((x*z*dx+y*z*dy)/r**2/sqrt(x**2+y**2)-sqrt(x**2+y**2)*dz)/r**2\n",
" dph = -y/(x**2+y**2)*dx+x/(x**2+y**2)*dy\n",
" return (dr, dth, dph)"
]
},
{
"cell_type": "code",
"execution_count": 41,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[ -cos(beta)*cos(theta)*sin(alpha) - (sin(alpha)*sin(beta)*sin(phi) - cos(alpha)*cos(phi))*sin(theta) r*cos(beta)*sin(alpha)*sin(theta) - (sin(alpha)*sin(beta)*sin(phi) - cos(alpha)*cos(phi))*r*cos(theta) -(cos(phi)*sin(alpha)*sin(beta) + cos(alpha)*sin(phi))*r*sin(theta)]\n",
"[ cos(beta)*sin(phi)*sin(theta) - cos(theta)*sin(beta) r*cos(beta)*cos(theta)*sin(phi) + r*sin(beta)*sin(theta) r*cos(beta)*cos(phi)*sin(theta)]\n",
"[ cos(alpha)*cos(beta)*cos(theta) + (cos(alpha)*sin(beta)*sin(phi) + cos(phi)*sin(alpha))*sin(theta) -r*cos(alpha)*cos(beta)*sin(theta) + (cos(alpha)*sin(beta)*sin(phi) + cos(phi)*sin(alpha))*r*cos(theta) (cos(alpha)*cos(phi)*sin(beta) - sin(alpha)*sin(phi))*r*sin(theta)]\n"
]
}
],
"source": [
"def print_formulas(): # enclosed in a function to prevent altering the namespace\n",
" alpha, beta = var('alpha, beta')\n",
" spher. = E.chart()\n",
" x, y, z = r*sin(theta)*cos(phi), r*sin(theta)*sin(phi), r*cos(theta) #Â normal Spherical->Cartesian transformation\n",
" y, z = y*cos(beta)-z*sin(beta), z*cos(beta)+y*sin(beta) # first rotation\n",
" x, z = x*cos(alpha)-z*sin(alpha), z*cos(alpha)+x*sin(alpha) # second rotation\n",
" spher.transition_map(E.default_chart(), [x, y, z])\n",
" print(list(E.changes_of_frame().values())[0][spher.frame(),:, spher])\n",
"\n",
"print_formulas()"
]
},
{
"cell_type": "code",
"execution_count": 42,
"metadata": {},
"outputs": [],
"source": [
"%%cython\n",
"from libc.math cimport exp, erf\n",
"\n",
"cpdef float profile(float x, float disk_min, float disk_max):\n",
" cdef float y \n",
" # we really don't want negative values\n",
" if xdisk_max:\n",
" return 0\n",
" y = (exp(-(disk_min-20-x)**2/400)*(x-disk_min)**2*(disk_max-x)**2/10000 +\n",
" exp(-(32-x)**2/70)/2*(x-disk_min)**2*(disk_max-x)**2/150000)\n",
" return max(y, 0)"
]
},
{
"cell_type": "code",
"execution_count": 43,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAkwAAAGFCAYAAAAPa6wiAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAzlElEQVR4nO3deZzd0/3H8dfJYkIkg6SKmNilgjaW2AWtvYu1irZarX0NQSVtqbY6WpHYpvx+drXToopfxZYgMbFFEVuLCmJNzUTIkOT7++NM0jEyubPcO+cur+fj8X3czPfeO99PjpF533PO95yQZRmSJElqW4/UBUiSJBU7A5MkSVIOBiZJkqQcDEySJEk5GJgkSZJyMDBJkiTlYGCSJEnKwcAkSZKUQ0kEphD1DyGE1LVIkqTK0yvx9du1zHhDQwPV1dU0NDQUuh5JklRZ2tUZUxI9TJIkSSl1ODCFEEaEEO4IIbwVQshCCHt24L1bhxDmhRCmdfS6kiRJqXSmh6kv8DRwTEfeFEKoBq4G7uvENSVJkpLp8BymLMvuBu4G6OAc7P8BrgPmA3t29LqSJEmpdMscphDCwcBawBndcT1JkqR8KvhdciGEdYCzgG2zLJvXnl6ppqYmmpqaFn3d2NhYuAIlSZJyKGgPUwihJ3EY7vQsy15q7/tqa2uprq5edNTU1BSuSEmSpBxClrVrKaTFvzmEDNgry7Lb2nh+OeA/xHlLC/Ugrnkw/7777uv59a9//QvvW1wPU01NDQ0NDfTv37/T9UqSJLXSrgnZhR6SawQ2bHXuKODrwL6bb775M4t7U1VVFVVVVQUuTZIkqX06HJhCCMsCa7c4tUYIYRgwK8uy10MItcCgLMsOyrJsAfBsq/e/C8zNsuxz5yVJkopVZ3qYNgUeaPH1uObHq4AfAysDg7tWliRJUvHo0hymPGjXxRsbGxftJeccJkmSlEfuJSdJkpQPBiZJkqQcDEySJEk5GJgkSZJyMDBJkiTlYGCSJEnKoagDU11dHUOHDmX48OGpS5EkSRXMdZgkSVIlcx0mSZKkfDAwSZIk5WBgkiRJysHAJEmSlIOBSZIkKQcDkyRJUg4GJkmSpBwMTJIkSTkYmCRJknIwMEmSJOVQ1IHJveQkSVIxcC85SZJUydxLTpIkKR8MTJIkSTkYmCRJknIwMEmSJOVgYJIkScrBwCRJkpSDgUmSJCkHA5MkSVIOBiZJkqQcDEySJEk5FHVgci85SZJUDNxLTpIkVTL3kpMkScoHA5MkSVIOBiZJkqQcDEySJEk5GJgkSZJyMDBJkiTl0OHAFEIYEUK4I4TwVgghCyHsmeP1e4cQJoQQ3gshNIYQpoQQdul0xZIkSd2sMz1MfYGngWPa+foRwARgd2AT4AHgjhDCRp24tiRJUrfr1dE3ZFl2N3A3QAi513rKsmxkq1NjQgh7AN/u6LUlSZJS6HBg6qoQQg+gHzCrrdc0NTXR1NS06OvGxsZuqEySJGnxUkz6HkUc1ruprRfU1tZSXV296Kipqem+6iRJklrp0l5yIYQM2CvLstva+foDgEuBPbIsu5c29pJbXA9TTU2Ne8lJkqR8a9dect02JBdC+B5wGfDd5rDUpqqqKqqqqrqnMEmSpBy6ZUiuuWfpSuDALMvu7I5rSpIk5UuHe5hCCMsCa7c4tUYIYRgwK8uy10MItcCgLMsOan79AcDVwPHAoyGElZrf90lXhgMlSZK6S2eG5DYlrqW00Ljmx6uAHwMrA4NbPH9483Xqmg9avF6SJKnodWnSdx606+KNjY1UV1c76VuSJOVbuyZ9u5ecJElSDgYmSZKkHAxMkiRJORiYJEmScjAwSZIk5VDUgamuro6hQ4cyfPjw1KVIkqQK5rICkiSpkrmsgCRJUj4YmCRJknIwMEmSJOVgYJIkScrBwCRJkpSDgUmSJCkHA5MkSVIOBiZJkqQcDEySJEk5GJgkSZJyKOrA5F5ykiSpGLiXnCRJqmTuJSdJkpQPBiZJkqQcDEySJEk5GJgkSZJyMDBJkiTlYGCSJEnKwcAkSZKUg4FJkiQpBwOTJElSDgYmSZKkHIo6MLmXnCRJKgbuJSdJkiqZe8lJkiTlg4FJkiQpBwOTJElSDgYmSZKkHHqlLkAqNvPnQ3093HUX3H03zJgBn3wSn1t3XVh//XjsuCNsvDGEdk0XlCSVMu+Sk5q99BKceSb87W8waxYMGAC77BLD0dJLw7x58MILMH06PPssfPQRrLEGHHEEHHIIrLBC6r+BJKkT2vWxt8OBKYQwAjgZ2ARYGdgry7LbcrxnO2AcsD7wFvCHLMsuxsCkIvDRRzB6NFx0EQwaBD/6Eey+OwwfDj17Lv49n30GEyfCNdfA9ddDr15wyilw8smwzDLdW78kqUsKtqxAX+Bp4Jh2VRHCGsBdwEPARsDvgPNDCPt04tpSXk2fDpttBldcAbW18OKL8OtfwxZbtB2WAHr3jkNyV14Zh+yOOQZ+9ztYbz246SZI23ErScq3Lg3JhRAycvQwhRB+D3wny7L1Wpy7GPhalmVbtOc69jCpEG64AX76U1h9dbjllhh2uuKf/4RRo+Cvf4Xtt4erroLBg/NRqSSpgIpm4cotgXtanfs7sOlnn3222Dc0NTXR2Nj4uUPKp8svhwMOgD33hKlTux6WANZeG26/Hf7+d/jXv2DYMLj11q5/X0lSet0RmFYC3ml17h2g1/vvv7/YN9TW1lJdXb3oqKmpKXSNqiBXXhknaR95ZJyD1Ldvfr//zjvDtGmxl2nvveNw3dy5+b2GJKl7ddc6TK3H/QJAaON+7NGjR9PQ0LDomDFjRqHrU4X485/hJz+BQw+FCy8s3JIAK6wQr1VXB5deCjvsAO++W5hrSZIKrzsC09vEXqaWVgTmDRgwYLFvqKqqon///p87pK569tl4B9x3vxvviOtR4J/+EOCoo+Chh+C11+JE8unTC3tNSVJhdEdgmgLs1OrczsDjvXv37obLS/Dhh7DXXrDmmnH+UqHDUkvDh8eFMPv2ha22gnvv7b5rS5Lyo8O/NkIIy4YQhoUQhjWfWqP568HNz9eGEK5u8ZaLgdVCCONCCOuFEH4C/BQY29XipfbIMvjhD+H99+Mk7HzPWWqPwYPhkUdiL9Nuu8G113Z/DZKkzuvM5+xNgaeaD4gLUj4F/Lr565WBRTdTZ1n2KrA7sD0wDfglcFyWZX/uVMVSB112WVy9+5prYK210tXRv3+s44c/hB/8IM5vkiSVBrdGUVl7800YOhT23TcGp2KQZXDSSTBuXNyKZcyY1BVJUkVr1+0/br6rspVlcZ+3vn3hnHNSV/NfIcDYsVBdDT//eTxnaJKk4mZgUtm6/vo4BHb77bDccqmr+bwQ4LTT4p9//vP49ejRaWuSJLXNwKSyNGdO3KZkv/3gO99JXU3bTjst9oSNGRND06mnpq5IkrQ4BiaVpfHjYdYs+P3vU1eS2+mnx9A0enRc7uCUU1JXJElqzcCksvPee/CHP8QtSVZfPXU17XP66bBgAfzsZ7Gn6eSTU1ckSWqpqANTXV0ddXV1zJ8/P3UpKiFnnhlDRylNpA4Bzjgj9jSdckr8+qSTUlclSVrIZQVUVl59FYYMieGjFCdRZxn88pcx9I0dG+dhSZIKymUFVHl++UsYOBCOPz51JZ0TAvzmN/9dq6lHDzjhhNRVSZIMTCobzzwTtxz53/+FZZZJXU3nhQC//W2c03TiidCrFxx7bOqqJKmyGZhUNs45J+7ZdvDBqSvpuhDgd7+Dzz6D446D3r3jIpySpDQMTCoLb78dF6o888zYI1MOQoCzz4Z58+DII+Pf65BDUlclSZWpTH61qNJddFHshSm3QBFCXFPqs8/gsMPi3/FHP0pdlSRVHgOTSt7cuTEwHXxw8W2Bkg8hwAUXxNB08MGxp+n7309dlSRVFgOTSt6118L778e5PuWqRw+4+OI4PHfQQTE0fe97qauSpMphYFJJy7I4ZPWtb8E666SuprB69IBLLomhaWEPk6FJkrqHgUkl7b774Lnn4pBVJejZE664Iv75wAPhk0/gxz9OWpIkVQQDk0ra+PHwta/B9tunrqT79OwJV14JSy8d5zR9/DEcdVTqqiSpvBmYVLJeeAHuuiv2uIR2LWxfPhbOaVpmGTj66Bia3HtOkgqnqAOTm+9qSS68EL78ZTjggNSVpBECjBsHffvCySfDnDlw2mmVFx4lqTu4+a5K0ty5sPLKcfXr2trU1aRXWwtjxsRept//PvZASZLaxc13Vb7++lf48EMnPC80ejQsu2zcdHjGjDjHqU+f1FVJUvnwc6hK0pVXwpZbwpAhqSspHsceCzffDLffDjvuGNemkiTlh4FJJeett+Dvf3eLkMXZZx944AF46aUYKF9+OXVFklQeDEwqOddeG/dUc9HGxdtiC3j00bj8wJZbwkMPpa5IkkqfgUklJcvicNxee5XnvnH5suaaMHkybLABfP3rUFcX206S1DkGJpWUxx+H6dOd7N0eK6wAEybERS2POSbuQTd7duqqJKk0GZhUUq68ElZZJU5qVm69e8N558E118Ctt8Imm8CTT6auSpJKj4FJJeOzz+DGG+EHP4jzc9R+3/8+PPUU9OsX5zidfXbcxFeS1D4GJpWM++6DDz6o3JW9u2qddeK8puOOg5/9LE4If/rp1FVJUmkwMKlk3HADrLtu3GxXnVNVBWPHwpQpcbX0TTeFn/88/lmS1LaiDkx1dXUMHTqU4cOHpy5FiTU1xTk4++/vXmn5sPnm8MQTce+5sWNhww3hz3/2TjpJaot7yakk/PWvsMce8NxzMHRo6mrKy/PPw6hRcPfd/53ftM02qauSpG7Tro/hRd3DJC10662w3nqGpUJYbz246y64997Yk7fttrDrrjBpUurKJKl4GJhU9ObNgzvugD33TF1JefvGN+I6V9dfD2++CdttFyeGX3cdfPpp6uokKS0Dk4reQw/Fu+P22it1JeWvR484T+zpp+MwaN++cUmCmho49dS4aKgkVSIDk4rebbfBoEHxji51jx494NvfjsN006fDfvvB//4vrL9+XPyytjbOfXKSuKRK4aRvFbUsg9VXj7+8L7wwdTWVrakpznW69lr4v/+DOXPi2k577AE77QRbbx17pCSpxBRu0ncI4agQwqshhLkhhCdCCNvmeP33QwhPhxA+DiHMDCFcEUIY0Jlrq7I89RS8/rrDccWgqir+d7jlFnjvvTivbLvt4rYru+wCyy8PI0bAr34FEyfGgCVJ5aLDPUwhhO8BfwKOAh4BDgcOAYZmWfb6Yl6/DTAROAG4AxgEXAy8nGXZnu25pj1MleuXv4S6OnjnnbgvmopPlsELL8D998fjgQfgP/+JAWuLLWKoGjEiTiBfZpnU1UrSF7Srh6kzgakeeDLLsiNbnHseuC3LstGLef1JwJFZlq3V4tyxwClZlq3anmsamCrXBhvAxhvD1VenrkTttWBBnDQ+cWI8Jk2CWbNi4B0+PIan7baLQ3j9+qWuVpIKEJhCCEsBHwPfzbLs1hbnzwOGZVm23WLesxXwALAXcDewInAT8HyWZYe357oGpsr08stxK5S//MUhuVK2YEGcON4yQL3zTpxYvvHGcf7TN78Ze6B6eBuKpO5XkDlMA4GewDutzr8DrLS4N2RZNhn4PnAj8CnwNvAhcGxbF2lqaqKxsfFzhyrPbbfB0kvH+TEqXT16xJ7Co4+Gm26CmTPjEN5FF8VJ45dcElcWX3VVOPbYGKgWLEhdtSR9Xmc/z7XulgqLORefCGEocD7wa2ATYFdgDeI8psWqra2lurp60VFTU9PJMlXKbr8ddt7ZeS/lJgQYMgQOOywuivn223Gtrf32iyu6b7cdrL02nHlmDFeSVAy6Y0juT0CfLMu+2+LcNsBDb731FiuvvPIXrtPU1ERTi1tsGhsbqampcUiugsyaBV/6EvzP/8Ahh6SuRt1lwQJ45BG4/PLYGzVvHvz0p/Czn8Fqq6WuTlKZyv+QXJZlnwJPADu1emonYHIbb1sGaN3BPr/5+y32DVVVVfTv3/9zhyrLPffEX5677Za6EnWnHj3iXnZXXAFvvRWXKLj55tjjdMgh8M9/pq5QUqXqzJDcOOCQEMJPQgjrhRDGA4NpHmILIdSGEFre03QHsHcI4cgQwpohhK2JQ3RTV1llla7WrzJ1113wta/FFb5VmaqrYfRoeO01OOss+Nvf/juU9957qauTVGk6HJiyLLsRGAmcBkwDRgC7Z1n27+aXrEwMUAtffyVwInAM8CxwM/AisHfny1Y5W7AgriS9++6pK1Ex6NsXRo2CV1+FceNij9O668JVV7k1i6Tu49YoKjqPPQabbRYnAm+zTepqVGzefx9OPBH+9CfYZ584z22A+wZI6rzCbY0iFdJdd8Fyy8VVoqXWBg6MC5nefHNcVXzDDeHvf09dlaRyZ2BS0bnrrrj2Uq9eqStRMdt3X3jmmRiYdt0VTjoJ5s9PXZWkcmVgUlF57704JOfdcWqPVVaBu++Oc5vGj48h6uOPU1clqRwZmFRU7rsvTuTdeefUlahU9OgBJ5wQFzq95x7Yfvu49Yok5ZOBSUVlwoS4jcZi1jOVluhb34o3CrzxBmy+edy/TpLyxcCkopFlMTDtuGPqSlSqNt4Y6uuhX7+4xcqzz6auSFK5KOrAVFdXx9ChQxk+fHjqUtQNXn4ZZsyIu9dLnVVTAw8+GDfz3WEHQ5Ok/HAdJhWNuro4F2XWLFh22dTVqNR98EHsrXzzTXj44bjYpSQthuswqbTcey9suaVhSfkxYED8mRo4MN5E8NZbqSuSVMoMTCoK8+bB/fc7HKf8GjAgLmo5f35cq6mxMXVFkkqVgUlF4bHH4i8zJ3wr32pqYmh6/XU48EAXt5TUOQYmFYV7742702+6aepKVI6GDoUbb4yLXI4Zk7oaSaXIwKSiMGECfP3rboeiwtllFxg7Fv7wh7hxryR1hIFJyc2eDVOmOBynwhs5Eg4+GA49FKZNS12NpFJiYFJykybFSd9O+FahhQB//COstx7svz/MmZO6IkmlwsCk5B58EAYNgrXXTl2JKkGfPnDDDXGR1OOOS12NpFJhYFJyEyfGbSxCu5YOk7puyBC48EK4/PIYniQpFwOTkpo9G558MgYmqTv9+MdwwAFw+OHwyiupq5FU7Io6MLmXXPl75JG4Lo6BSd0tBLjoori45QEHwGefpa5IUjFzLzklNXo0XHEFzJzpkJzSqK+HrbaC3/zGNZqkCuVecip+kyY5f0lpbb45nHIKnHEGPPdc6mokFSsDk5L5+OO4JYrDcUrt9NNhrbXiGk3z5qWuRlIxMjApmSlT4ryRESNSV6JK16dPHBp+/HG44ILU1UgqRgYmJTNxYpxwO3Ro6kqkODR31FFw2mnwxhupq5FUbAxMSmbixNi71MOfQhWJ3/4W+vaFE05IXYmkYuOvKiUxd268O8n5Syomyy0H48bBLbfA//1f6mokFRMDk5Kor4emJgOTis8BB8A3vgFHHw2ffJK6GknFwsCkJCZNip/mN9wwdSXS54UAdXVxHlNtbepqJBULA5OSmDgRttkGevZMXYn0RUOGwKhRcPbZ8PrrqauRVAwMTOp2n34Kkyc7HKfiNnp07AU99dTUlUgqBkUdmNxLrjw9/nicG2JgUjHr1w/OPBOuvz6uGSapsrmXnLpdbW08Zs2CXr1SVyO1bf582HRTqKqKvaIugSGVJfeSU3GaOBG23tqwpOLXsyeMHx/v6rzhhtTVSErJwKRuNW8ePPKIw3EqHdtvD3vsAb/4RZx/J6kyGZjUrZ5+Gj76CLbdNnUlUvv99rfw2mtw6aWpK5GUioFJ3WryZFhqKdhkk9SVSO23wQbwgx/Ab34Dc+akrkZSCgYmdavJk2NY6tMndSVSx5xxBnzwAVxwQepKJKVgYFK3mjwZttoqdRVSx62xBhx+OPz+9/Cf/6SuRlJ361RgCiEcFUJ4NYQwN4TwRAhhiTNSQghVIYQzQwj/DiE0hRD+FUL4SedKVql64424avLWW6euROqchRO///CH1JVI6m4dDkwhhO8B5wJnAhsBDwF3hxAGL+FtNwHfAH4KDAEOAF7o6LVV2iZPjo9bbpm2DqmzvvxlGDkSzjsPZs5MXY2k7tSZHqYTgcuyLLs0y7LnsywbCcwAjlzci0MIuwLbAbtnWXZvlmWvZVk2NcuyyZ2uWiVp8mRYc01YaaXUlUidd/LJcQ7eb36TuhJJ3alDgSmEsBSwCXBPq6fuAdqamfId4HHglBDCmyGEl0IIY0MIS7d1naamJhobGz93qPQ5f0nlYLnl4JRT4LLL4jCzpMrQ0R6mgUBP4J1W598B2uo3WBPYBtgA2AsYCewL1LV1kdraWqqrqxcdNTU1HSxTxebjj+Gpp5y/pPJw9NHQt69zmaRK0tm75FrvARcWc67lNTLg+81DcXcRh/V+/Mknnyz2DaNHj6ahoWHRMWPGjE6WqWLx2GNxlW97mFQO+vWDE06ASy6Bt99OXY2k7tDRwPQ+MJ8v9iatyBd7nRaaCbyZZVlDi3PPA+GNNvqzq6qq6N+//+cOlbbJk+MvmfXXT12JlB/HHhsXYR07NnUlkrpDhwJTlmWfAk8AO7V6aiegrUncjwCrhBCWbXFuXWDBqquu2pHLq4RNngxbbBE3M5XKwXLLwXHHwUUXwXvvpa5GUqF1ZkhuHHBICOEnIYT1QgjjgcHAxQAhhNoQwtUtXn8d8AFwRQhhaAhhBHA2cPnSS7c571tlJMtiYHL+ksrNyJHQoweMG5e6EkmF1uHAlGXZjcSJ26cB04ARxCUD/t38kpWJAWrh6z8i9kAtR7xb7lrgDuC4zpetUvLiizBrlvOXVH4GDICjjoILL4w/45LKV8iytuZqd4t2XbyxsZHq6moaGhqcz1SCLr8cDjkEPvwQ/M+ncvPuu7D66nDSSfDrX6euRlInhPa8yL3kVHCTJ8OGGxqWVJ5WXBGOOALOPx8aGnK/XlJpMjCp4FywUuXu5JNh7ly44ILUlUgqFAOTCmrWLHj+eSd8q7ytvHIcdj733LhIq6TyY2BSQU2ZEh/tYVK5GzUK/vMfuPLK1JVIKgQDkwpqypQ4x2ONNVJXIhXWGmvAvvvCOefA/Pmpq5GUb0UdmOrq6hg6dCjDhw9PXYo6aepU2HxzCO26B0EqbSefDK+8An/5S+pKJOWbywqoYBYsgBVWiL9Efv7z1NVI3WOHHWDOHKiv94OCVCJcVkBpvfxyvM16s81SVyJ1n5NPjptNT5qUuhJJ+WRgUsHU18dHR1RVSXbbDTbYAM4+O3UlkvLJwKSCqa+Hr3wlblIqVYoQ4qrfd94J06enrkZSvhiYVDBTpzocp8p0wAEwaBCMHZu6Ekn5YmBSQcydC08/He+QkyrNUkvB8cfDNdfAW2+lrkZSPhiYVBDTpsFnn9nDpMp12GHQp0/cY05S6TMwqSDq66GqCr761dSVSGlUV8Phh8PFF0NjY+pqJHWVgUkFMXUqbLRRHJqQKtXxx8c1mS65JHUlkrrKwKSCqK93/pK06qpw4IFxU97PPktdjaSuMDAp795/H/71LwOTBHGJgTfegBtvTF2JpK4o6sDkXnKl6bHH4qMTviXYcEPYaSc47zxIuxOVpK5wLznl3a9+BRdeCO+9515aEsDdd8Puu8PDD8PWW6euRlIr7iWnNBYuWGlYkqJddoEhQ2D8+NSVSOosA5PyKstiYHL+kvRfPXrEO+ZuvRVeey11NZI6w8CkvHrlFfjgA+cvSa0ddBD07x+HqyWVHgOT8qq+Pj4amKTP69s3rv596aUwe3bqaiR1lIFJeTV1Kqy9NgwYkLoSqfgccwx89BFcdVXqSiR1lIFJeVVfb++S1JaaGth337jEwIIFqauR1BEGJuXNp5/CU0854VtakpEj4Z//hDvvTF2JpI4wMClv/vEPaGqyh0laki22iB8qzj03dSWSOsLApLypr4fevWHYsNSVSMVt5Ei4//74IUNSaTAwKW+mToWvfQ369EldiVTc9tknbsx73nmpK5HUXkUdmNxLrrTU1zt/SWqP3r3jHXPXXgvvvpu6Gknt4V5yyosPP4Tll4err4Yf/jB1NVLxmzUr9jKdeiqcdlrqaqSK5l5y6j6PPRYfnfAttc8KK8CPfgR//GO8WUJScTMwKS/q62G55WCddVJXIpWO44+Hd96BG29MXYmkXAxMyoupU2H48LjJqKT2+cpXYLfdYPz4uHG1pOLlrzd1WZY54VvqrJEjYdo0mDgxdSWSlsTApC57/fV4p4/zl6SO22knWG89lxiQip2BSV1WXx8f7WGSOi6E2Mt0++3wyiupq5HUlk4FphDCUSGEV0MIc0MIT4QQtm3n+7YOIcwLIUzrzHVVnKZOhdVXhxVXTF2JVJp+8IO4LMcFF6SuRFJbOhyYQgjfA84FzgQ2Ah4C7g4hDM7xvmrgauC+jpepYlZf73Cc1BXLLAOHHw6XXQaNjamrkbQ4nelhOhG4LMuyS7Msez7LspHADODIHO/7H+A6YEonrqkiNW8ePPGEw3FSVx19NHzyCVxxRepKJC1OhwJTCGEpYBPgnlZP3QNstYT3HQysBZzR0QJV3J59Nv4jbw+T1DWDBsF++8H558P8+amrkdRaR3uYBgI9gXdanX8HWGlxbwghrAOcBXw/y7J57blIU1MTjY2NnztUnOrroWdP2Hjj1JVIpe/44+PE77/9LXUlklrr7F1yrZdYC4s5RwihJ3EY7vQsy15q7zevra2lurp60VFTU9PJMlVoU6fChhvGORiSumazzWCrreDcc1NXIqm1jgam94H5fLE3aUW+2OsE0A/YFLiw+e64ecBpwNdCCPPuv//+xV5k9OjRNDQ0LDpmzJjRwTLVXVywUsqvkSPhwQfjYpaSikeHAlOWZZ8CTwA7tXpqJ2DyYt7SCGwIDGtxXAy8CAzbvI3ftFVVVfTv3/9zh4rP7NkwfbqBScqnvfaCmhoXspSKTWeG5MYBh4QQfhJCWC+EMB4YTAxChBBqQwhXA2RZtiDLsmdbHsC7wNwsy57t27dvvv4eSuDxx+O2KE74lvKnVy849li47rq4Ma+k4tDhwJRl2Y3ASOLQ2jRgBLB7lmX/bn7JysQApTI3dSr06xc3EJWUP4ccEoPTRRelrkTSQiFLu0V2uy7e2NhIdXU1DQ0NDs8Vkb33hg8/hDamoknqgqOPhltugX//G/r0SV2NVNZCe17kXnLqNCd8S4Vz3HFxU+sbbkhdiSQwMKmT3nwT3nrL+UtSoQwZAt/8JowfH+cKSkrLwKROqa+Pj/YwSYVz4onwj3/AvfemrkSSgUmdMnUqrLoqrLJK6kqk8rXDDrDRRnDOOakrkWRgUqfU1zscJxVaCDBqFPz97/DMM6mrkSqbgUkdNn9+XIPJ4Tip8PbbL/bmjhuXuhKpshV1YKqrq2Po0KEMHz48dSlq4fnn4aOP7GGSukPv3nFT3muvjTdaSErDdZjUYZddBocdBg0NsOyyqauRyl9DQ9wu5eijobY2dTVS2XEdJhXG1KkwdKhhSeou1dVw6KFw8cWxd1dS9zMwqcNcsFLqfscfHze8vvzy1JVIlcnApA6ZMweefdbAJHW3wYPjBPBzz4V581JXI1UeA5M65Mkn411yTviWut+oUfDqq3DrrakrkSqPgUkdMnUqLLMMrL9+6kqkyrPJJnExy7Fj3S5F6m4GJnVIfX38R7tXr9SVSJVp1Kj4weXhh1NXIlUWA5M65NFHYYstUlchVa7ddot3qf7+96krkSqLgUnt9uabMGOGgUlKqUcPOPVUuPNOePrp1NVIlcPApHarr4+PBiYprf33h9VWg7POSl2JVDkMTGq3Rx+Nqw2vskrqSqTK1rs3nHIK3HQT/POfqauRKkNRByb3kisu9fX2LknF4uCDYeBAOPvs1JVIlcG95NQu8+ZB//7w29/CiSemrkYSxCG500+PazPZ8yt1mnvJKX+eeQY++cQVvqVicuSRsPTSMG5c6kqk8mdgUrs8+mhce2njjVNXImmh6mo4+ui4Ke8HH6SuRipvBia1y6OPwrBh8dOspOJx/PGwYAFceGHqSqTyZmBSu7hgpVScVlwRDjkEzj8fPvoodTVS+TIwKadZs+CllwxMUrE66SSYPdteJqmQDEzKaerU+OiEb6k4DR4Mhx4alxhobExdjVSeDEzK6dFHYcAAWGut1JVIasuYMTBnDpx3XupKpPJkYFJOC+cvhXatVCEphUGD4Igj4Jxz4MMPU1cjlR8Dk5ZowQJX+JZKxamnwqefui6TVAgGJi3RSy/FT6sGJqn4rbRSXJfp3HNdl0nKt6IOTO4ll96jj8ahOP8TSKXhlFNiz/DYsakrkcqLe8lpiY44Ah56CJ57LnUlktprzJi4LtMrr8R1miQtkXvJqeucvySVnpNOgh494jIDkvLDwKQ2zZ4N//iHgUkqNSusACecAHV18PbbqauRyoOBSW2qr49zIbbeOnUlkjrqhBOgqgrOOit1JVJ5MDCpTZMnw/LLw1e+kroSSR213HIwahRcfDG88UbqaqTSZ2BSmx55BLbaKs6FkFR6jj8e+veHX/4ydSVS6fNXoRZr/nyYMiUGJkmlqV8/OOMMuOoqmDYtdTVSaetUYAohHBVCeDWEMDeE8EQIYdslvHbvEMKEEMJ7IYTGEMKUEMIunS9Z3eG55+Kkb+cvSaXt0ENhyJA4PJd2FRmptHU4MIUQvgecC5wJbAQ8BNwdQhjcxltGABOA3YFNgAeAO0IIG3WmYHWPRx6BXr1csFIqdb16xeUF7r8f7rwzdTVS6erwwpUhhHrgySzLjmxx7nngtizLRrfzezwH3Jhl2Rnteb0LV3a/H/4QXnwRpk5NXYmkrsoy2HFHeOutuFRI796pK5KKSv4XrgwhLEXsJbqn1VP3AO2a7RJC6AH0A2a19ZqmpiYaGxs/d6h7PfKIw3FSuQghbpXy4otw6aWpq5FKU0eH5AYCPYF3Wp1/B1ipnd9jFNAXuKmtF9TW1lJdXb3oqKmp6WCZ6oqZM+HVV53wLZWTjTaCgw6C008HP4NKHdfZu+Raj+OFxZz7ghDCAcCvgO9lWfZuW68bPXo0DQ0Ni44ZM2Z0skx1xqRJ8XHbNqfySypFv/0tfPQRnHlm6kqk0tPRwPQ+MJ8v9iatyBd7nT6nebL4ZcB+WZbdu6TXVlVV0b9//88d6j6TJsG668JK7e0zlFQSVl0VTj0Vxo2DZ59NXY1UWjoUmLIs+xR4Atip1VM7AZPbel9zz9KVwIFZlnmfRpGbNAlGjEhdhaRC+NnPYK214Igj4tZHktqnM0Ny44BDQgg/CSGsF0IYDwwGLgYIIdSGEK5e+OLmsHQ1ce7SoyGElZqP6jzUrzz74IP4ydPAJJWnqiq46KJ4Y8dll6WuRiodHQ5MWZbdCIwETgOmEddZ2j3Lsn83v2RlYoBa6HCgF1AHzGxxnNfZolU4Dz8cHw1MUvnaYQf40Y9ib9O7bc4mldRSh9dhyrN2Xdx1mLrPqFHw5z/Da6+lrkRSIb3/flwBfPfd4U9/Sl2NlFT+12FS+Zs40d4lqRIMHBjXZrrmGrh3ibfhSAIDk1pobISnnjIwSZXixz+O/78fdRTMnZu6Gqm4GZi0yOTJ8a4ZA5NUGUKAiy+OQ/C1tamrkYqbgUmLTJoEX/4yrLNO6kokdZf11ouTv2tr4YUXUlcjFS8DkxZ54AHYbrv4qVNS5RgzBgYPhsMOg/nzU1cjFaeiDkx1dXUMHTqU4cOHpy6l7H34IUydCju1XpJUUtlbeum4JtPDD8P48amrkYqTywoIgL/8BfbZJ85lWG211NVISmHUKLjwQnjiCdhgg9TVSN3GZQXUfhMmxLlLhiWpcp15Zvx34Ac/gKam1NVIxcXAJCAGJofjpMrWp09cl+n55+HEE1NXIxUXA5N45RX4179g551TVyIptWHD4Pzz4Y9/hOuuS12NVDwMTGLCBOjZE7bfPnUlkorBYYfFYbnDDoPp01NXIxUHA5OYMAE23xyqq1NXIqkYLFzQcvXVYd994aOPUlckpWdgqnDz58N99zl/SdLn9e0Lt9wCr78ee5rS3lAtpWdgqnCPPx7XYHL+kqTWvvIVuPRSuP762OMkVTIDU4WbMAH694fNNktdiaRitP/+cMwxMHIk1NenrkZKx8BU4e65B3bYAXr1Sl2JpGI1dixsuil85ztxcVupEhmYKtjs2TBlisNxkpasqgpuuw2WXRa++c04jC9VmqIOTO4lV1gTJ8K8eU74lpTbl74Ed90FM2fGbZRcCVyVxr3kKtjhh8c75F5+Od5GLEm5TJoEu+wCu+8ON97ocL7KgnvJqW0LFsBf/wp77GFYktR+I0bAzTfD7bfH5QYWLEhdkdQ9DEwV6rHH4O23Y2CSpI741rfgqqvgyivhpz+NQ/tSubMztULdfjsMGABbbZW6Ekml6Pvfj73TBx0UVwK/9lpYaqnUVUmFYw9Thbr99vgp0fkHkjrrwAPhL3+BO+6IvdUff5y6IqlwDEwV6IUX4oaae+6ZuhJJpe4734E774SHHoJdd4WGhtQVSYVhYKpA118fV/feddfUlUgqB9/4Rtw14Jln4p/ffTd1RVL+GZgqTJbBDTfA3ntDnz6pq5FULrbcEh58EN54A4YPhyefTF2RlF8Gpgrz1FPw0ktxfyhJyqevfS3egbviirD11nDddakrkvLHwFRhbrgBBg6M3eaSlG81NXFxy/32i3fSHXkkzJmTuiqp6wxMFeSTT+LaKfvv791xkgpn6aXjGk0XXRT/zdlkE3jiidRVSV1T1IHJveTy6/LL4f334YQTUlciqdyFAEccEecy9e0LW2wBZ54Jn36aujKpc9xLrkJ89hmsvTZss01cYE6Susunn8Lpp8Mf/gDrrgsXXAA77pi6KmkR95LTf11/Pbz+Opx6aupKJFWapZaC2tp408mXvgQ77QTf/S7MmJG6Mqn9DEwVYMECOOss+Pa3YcMNU1cjqVJ99aswcSJccw08/HDsbRo5EmbOTF2ZlJuBqQLcfjs8/zyMHp26EkmVLoR499yLL8Ye76uugjXWgOOOi2s4ScXKOUxlLstg881hmWXionKSVEwaGuKcpnHjoLEx7nF56KFxJ4KePVNXpwrhHCbB/ffHheTGjEldiSR9UXU1/OIX8NprMTi9/noMTautBqedFs9LxcAepjL3jW/Ahx/C44/HrnBJKnZPPAGXXBJXCp89G4YNi3Mwv/Ut2HRT6OFHfeVXu347GpjK2OTJcXuCm2+GffdNXY0kdcycOXDHHfC3v8Fdd8F//gNf/jLssgtsu2081l3XD4PqMgNTJfvgg/hJbOBAePRR5wJIKm3z5sGUKTFA3XcfTJsW7wDu3z/uYTds2H8fhwyBZZdNXLBKiYGpUr34YrwL5d//jl3bgwenrkiS8quxMX4YfPLJGJ6efjr+27fwV9rAgbD66vFYbbX4OGgQDBgQj4EDYYUVoHfvdH8HFY12BaZkO4qFEEJDQ8Nin2tqaqKpqWnR19dd9wkA55/fSJ8+3VJeyViwIK6i29QUj4aGeJvuoEFwyy2w3HLxHxZJKjdbbBGPhT7+GKZPh3/9K04eX3g8+WR8nDfvi9+jXz9Yfvm4fcsyy3zx6NMn9tAvPHr1ikfrcwsfu2N4sNDXKJchzo03jiMtuVRXV/cHZmc5epCS9TCFEPoDi09MkiRJ3ac6y7Ildi+kDEyhoaFhweKea93DNHPmTDbbbDOmT5/OoEGDOn3N4cOH89hjj3X6/fn8PsVUS2NjIzU1NcyYMaNLQ57F9Hcqllry1bb5qCVf36OYavFnt3Dfp9h+dvP1fYqllmJr32Jpl3x9j460b3V1dTXt6GFKNiSXq7DF6devX5d+sHr27JmXOVD5+D7FVMtC/fv3t30L9H262rb5qqWY2sWf3eKvBYrnZzdf36eYaoHiad9iapfu/rchV8/SQhW1msXRRx9dNN+nmGrJl2L6OxVTLflSbu1Sbm2br+9TTLXkSzH9nYqplnwpt3YpprZtqSTuknvjjTcWda2tuuqqha6p4ngXYuHYtoVl+xaObVtYtm9hdbB9y2drlKqqqs89Kr+qqqo4/fTTbd8CsG0Ly/YtHNu2sGzfwipE+5ZED5NJXJIkFUhJLFzZLi2WIMh5258kSVK+lUpgCkA/2nHbnyRJUr6VRGCSJElKqSQmfUuSJKVkYKoQIYQRIYQ7QghvhRCyEMKerZ4PIYRfNT//SQjhwRDC+onKLSkhhNEhhMdCCLNDCO+GEG4LIQxp9Rrbt5NCCEeGEP4RQmhsPqaEEHZr8bxtmyfNP8tZCOHcFuds305qbres1fF2i+dt2y4KIQwKIVwTQvgghPBxCGFaCGGTFs/nrY0NTJWjL/A0cEwbz58CnNj8/HDgbWBCCKFf95RX0rYD6oAtgJ2IK+jfE0Lo2+I1tm/nvQGcCmzafNwP3N7iHz3bNg9CCMOBw4B/tHrK9u2a54CVWxwbtnjOtu2CEMLywCPAZ8BuwFBgFPBhi5flr42zLPOosIO4nMOeLb4OwEzgZy3OVTX/0B2eut5SO4AvNbfxCNu3YG08C/ipbZu39lwWeAnYEXgQOLf5vO3btXb9FTCtjeds266371nAQ0t4Pq9tbA+TANYAVgLuWXgiy7ImYCKwVaqiSlh18+Os5kfbN09CCD1DCPsTe0ynYNvmSx1wZ5Zl97Y6b/t23TrNw0GvhhBuCCGs2Xzetu267wCPhxBubp4O8VQI4dAWz+e1jQ1MgvgDBfBOq/PvtHhO7dC8BMY44OEsy55tPm37dlEIYcMQwkdAE3AxsFeWZdOxbbusOYBuDIxezNO2b9fUAwcBuwCHEttscghhALZtPqwJHAm8TGzji4HzQwgHNT+f1zbu1ckiVZ5arzERFnNOS3Yh8FVgm8U8Z/t23ovAMGA5YB/gqhDCdi2et207IYRQA5wH7Jxl2dwlvNT27YQsy+5u8eUzIYQpwL+AHwGPLnxZq7fZtu3XA3g8y7IxzV8/1Ty38Ujg6havy0sb28MkiJPg4IuJe0W+mMzVhhDCBcQu4h2yLHujxVO2bxdlWfZplmX/zLLs8SzLRhNvYDge27arNiG21RMhhHkhhHnEmxiOa/7zwja0ffMgy7I5wDPAOvizmw8zgemtzj0PDG7+c17b2MAkgFeJP1g7LTwRQliK+A/n5FRFlYrm21YvBPYGvp5l2autXmL75l8gTt60bbvmPuJdW8NaHI8D1zb/+RVs37wJIVQB6xF/0fuz23WPAENanVsX+Hfzn/Paxg7JVYgQwrLA2i1OrRFCGAbMyrLs9eZ1V8aEEF4mjgePAT4GruvuWktQHXAgsAcwO4Sw8NNMQ5Zln2RZltm+nRdC+B1wNzCDuEXS/sD2wK62bddkWTYbeLbluRDCHOCDhXPwbN/OCyGMBe4AXif2avwC6A9c5c9uXownzgkbA9wEbEZcGuMwgHy3sYGpcmwKPNDi63HNj1cBPwb+ACwN/BFYnjhZcefmf1C1ZEc2Pz7Y6vzBwJXNf7Z9O+/LwJ+Ia9g0ENcJ2jXLsgnNz9u2hWX7dt6qwPXAQOA94rylLbIsW9gDYtt2QZZlj4UQ9gJqgdOIPUojsyy7tsXL8tbG7iUnSZKUg3OYJEmScjAwSZIk5WBgkiRJysHAJEmSlIOBSZIkKQcDkyRJUg4GJkmSpBwMTJIkSTkYmCRJknIwMEmSJOVgYJIkScrBwCRJkpTD/wONagNDas1XHwAAAABJRU5ErkJggg==\n",
"text/plain": [
"Graphics object consisting of 1 graphics primitive"
]
},
"execution_count": 43,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot(lambda x, d1=disk_min, d2=disk_max: profile(x, d1, d2), \n",
" xmin=0, xmax=60, ymin=0, ymax=1.4)"
]
},
{
"cell_type": "code",
"execution_count": 44,
"metadata": {},
"outputs": [],
"source": [
"%%cython\n",
"from libc.math cimport cos, sin, acos, sqrt, abs, atan2\n",
"cimport cython\n",
"from __main__ import profile\n",
"from __main__ import xyz_to_spherical\n",
"from __main__ import spherical_to_xyz\n",
"\n",
"@cython.boundscheck(False)\n",
"@cython.wraparound(False)\n",
"cpdef tuple intersection(list curve, float m, float alpha, float beta, \n",
" float dmin, float dmax):\n",
" cdef float x, y, z\n",
" cdef float x2, y2, z2\n",
" cdef float r, theta, phi\n",
" cdef int n, i\n",
" cdef float t\n",
" cdef float sinalpha, cosalpha\n",
" cdef float sinbeta, cosbeta\n",
" cdef float R, G, B\n",
" cdef float dr, dtheta, dphi\n",
" cdef float dx, dy, dz\n",
" cdef float th\n",
" R, G, B = 0., 0., 0. # return values\n",
" sinalpha = sin(alpha)\n",
" cosalpha = cos(alpha)\n",
" sinbeta = sin(beta)\n",
" cosbeta = cos(beta)\n",
" n = len(curve)\n",
" r, theta, phi = curve[0][2:5]\n",
" x, y, z = r*sin(theta)*cos(phi), r*sin(theta)*sin(phi), r*cos(theta)\n",
" x, y, z = x, y*cos(beta)-z*sin(beta), z*cos(beta)+y*sin(beta)\n",
" z = z*cos(alpha)+x*sin(alpha)\n",
" for i in range(1, n): \n",
" r = curve[i][2]\n",
" theta = curve[i][3]\n",
" phi = curve[i][4]\n",
" # rotations\n",
" x2, y2, z2 = r*sin(theta)*cos(phi), r*sin(theta)*sin(phi), r*cos(theta)\n",
" y2, z2 = y2*cosbeta-z2*sinbeta, z2*cosbeta+y2*sinbeta\n",
" x2, z2 = x2*cosalpha-z2*sinalpha, z2*cosalpha+x2*sinalpha\n",
" if z!=z2:\n",
" t = z/(z-z2)\n",
" if t>=0 and t<1 and curve[i][2]>dmin and curve[i][2]= 255:\n",
" R = 255\n",
" return R, G, B"
]
},
{
"cell_type": "code",
"execution_count": 45,
"metadata": {},
"outputs": [],
"source": [
"def render_row(x):\n",
" \"\"\"\n",
" Render a single row of the image\n",
" \"\"\"\n",
" res = np.zeros((ny,3)) # result row in RGB format\n",
" for y in range(ny):\n",
" beta = atan2(y-ny/2,x) # beta angle \n",
" r = sqrt(x**2+(y-ny/2)**2) # pixel distance to the center of the image\n",
" ind_geo = int(r/400*n_geod*720/nx) # index of the geodesic to use. values are obtained by trial and error.\n",
" if ind_geo"
]
},
"execution_count": 46,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"%display plain\n",
"data = np.zeros( (ny, nx, 3), dtype=float )\n",
"render()\n",
"img2 = toimage(data)\n",
"img2"
]
},
{
"cell_type": "code",
"execution_count": 47,
"metadata": {},
"outputs": [],
"source": [
"from six.moves.urllib.request import urlretrieve # valid for both Python 2 and Python 3\n",
"urlretrieve(\"http://www.cvrl.org/database/data/cmfs/ciexyzjv.csv\", \n",
" \"ciexyzjv.csv\")\n",
"ciexyz = np.genfromtxt(\"ciexyzjv.csv\", delimiter=\",\")"
]
},
{
"cell_type": "code",
"execution_count": 48,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAABIkElEQVR4nO3dd3yV1f3A8c+5N3uTvchmbwxTZCnLgaPu1boQUevP1tZVf7a1tlr9tbYuSi22rjoqBQcbGbIJyA6QCdl77+Se3x9PogEybpLnjtyc9+sVk9z7POccHvCbk+9ZQkqJoiiK0v8ZbN0ARVEURR8qoCuKojgIFdAVRVEchAroiqIoDkIFdEVRFAfhZKuKAwMDZUxMjK2qVxRF6ZcOHjxYLKUM6ug9mwX0mJgYkpKSbFW9oihKvySEONvZeyrloiiK4iBUQFcURXEQKqAriqI4CJvl0BVFUWyhqamJ7Oxs6uvrbd2ULrm5uREZGYmzs7PZ96iArijKgJKdnY23tzcxMTEIIWzdnA5JKSkpKSE7O5vY2Fiz71MpF0VRBpT6+noCAgLsNpgDCCEICAjo8W8RKqArijLg2HMwb9ObNnYb0IUQK4UQhUKI411cM1sIcVgIcUIIsb3HrXAw1Q3N/OdgNo3NJls3RVGUAcScHvo/gYWdvSmE8APeAhZLKUcBN+nSsn5KSsnPPz3ME58d4Z2d6bZujqIodui///0v48ePP+/DYDCwbt26PpXbbUCXUu4ASru45HZglZTyXOv1hX1qUT/37q5MNpwoIMTHlTe+SaWg0r5H0hVFsb7rr7+ew4cPf/+xbNkyLrvsMhYsWNCncvXIoQ8FBgkhtgkhDgoh7u7sQiHEEiFEkhAiqaioSIeq7cuhc2X8fm0y80aG8MmSaTS3SP64/rStm6Uoih07c+YMv/3tb3n//fcxGPoWkvWYtugEXAJcDrgDe4QQe6WUZy68UEq5AlgBkJiY6FBn35XXNvLoR98R6uvGqzeOw9fDmXtnxLJ8exp3TYtm/GA/WzdRUZQL/ObLE5zMrdS1zJHhPjx/zSizrm1qauL222/n1VdfJSoqqs9169FDzwbWSylrpJTFwA5gnA7l9iu/+zqZoqoG3rpjIr4e2kKAR+YmEOTtym++PIE6u1VRlAs999xzjBo1iltvvVWX8vTooa8B3hBCOAEuwBTgzzqU22+YTJJvThVy9bgwxkb6ff+6l6sTv1wwjF/85yhrDudy3YQI2zVSUZSLmNuTtoRt27bx+eefc+jQId3KNGfa4r+BPcAwIUS2EOI+IcRSIcRSACllMrAeOArsB96RUnY6xdERpRRWU1rTyLS4gIve+9HESIYEe/HR/nM2aJmiKPaorKyMe+65h/feew9vb2/dyu22hy6lvM2Ma14BXtGlRf3QvowSAKZ2ENANBsH8USEs355ORV0Tvu7m78ugKIpjWr58OYWFhTz00EPnvf70009zyy239LpctZeLDvamlxDh507kIPcO358zLJg3t6axM6WYq8aGWbl1iqLYm6effpqnn35a93LV0v8+klKyL72UKXH+nS7VHT/YD193Z745NaCn6CuKYmEqoPdRamE1JTWNTI29ON3SxsloYObQILafKcRkUrNdFEWxDBXQ+2hveuf58/bmDg+iuLqR47kV1miWoigDkArofbQ3vZQwXzcG+3ecP28zc0gQQsDWU463QlZRFPugAnofSCnZl1HC1Lju91YO8HJlXKQfW0+rPLqiKJahAnofpBVVU1zdyNQ4f7OunzMsmCPZ5ZRUN1i4ZYqiDEQqoPfBnnRtE8opXQyItjdneBBSwo4UlXZRlIFKSsmMGTPO2yr3008/ZeHCTncpN5sK6H2wN72EUB83ogM8zLp+dLgvgV6ufKPy6IoyYAkhWL58OT/72c+or6+npqaGZ599ljfffLPPZauFRb3UNv98RoL5ZxMaDILZw4LYdLKA5hYTTkb181RRBqLRo0dzzTXX8PLLL1NTU8Pdd99NfHx8n8tVAb2XMktqKa5uYLKZ6ZY2s4cF8Z+D2RzOKicxxrzcu6IoFrLuKcg/pm+ZoWNg0UvdXvb8888zceJEXFxcSEpK0qVqFdB7qW0P5TERvj2677KEIIwGwbbTRSqgK8oA5unpyS233IKXlxeurq66lKkCei8l51ViNAiGhHj16D5fD2cmRvmx7UwhTywYZqHWKYpiFjN60pZkMBj6fErReeXpVtIAk5xXSVygJ27Oxh7fO3tYMMdzKimsUueNKoqiHxXQeyk5r5IRYT69unf2sCAAtp9Ws10URdGPSrn0QnltI7kV9b0O6CPDfAj2dmXbmSJuShysc+sURekvfv3rX+tanjknFq0UQhQKIbo8hUgIMUkI0SKEuFG/5tmnU/lVAIwI691JI0IIZg0N4tszRTS3mPRsmqIoA5g5KZd/Al0uYRJCGIGXgQ06tMnuJedpM1xG9rKHDjBneDCV9c0czirXqVWKogx03QZ0KeUOoLSbyx4FPgcGxM5TyXmVBHi6EOTd+6lGlyYEYjQItVmXoii66fOgqBAiArgeWG7GtUuEEElCiKSiov47IJicV8WIMB+zV4h2xNfdmUuiBrFNDYwqiqITPWa5vAY8KaVs6e5CKeUKKWWilDIxKChIh6qtr7nFxOmCql7nz9ubNSyIE7lq+qKiKPrQI6AnAh8LITKBG4G3hBDX6VCuXcoorqGx2dTrGS7tzRkWDKB66Yqi6KLPAV1KGSuljJFSxgD/AZZJKVf3tVx7dbJ1QFSPgD4izJsIP3fWHcvrc1mKovQPWVlZxMbGUlqqDU2WlZURGxvL2bNn+1y2OdMW/w3sAYYJIbKFEPcJIZYKIZb2ufZ+KDmvCmejID6oZ0v+OyKE4OqxYXybUkx5baMOrVMUxd4NHjyYhx56iKeeegqAp556iiVLlhAdHd3nsrtdWCSlvM3cwqSUP+lTa/qB5LxK4oO8cHHSZ5Ht1WPD+duOdDacyOeWSVG6lKkoin17/PHHueSSS3jttdfYuXMnr7/+ui7lqpWiPZScV8mMhEDdyhsd4UN0gAdfHc1TAV1RrOzl/S9zqvSUrmUO9x/Ok5Of7PIaZ2dnXnnlFRYuXMjGjRtxcXHRpW61l0sPlFQ3UFjVoEv+vE1b2mV3Wok6a1RRBpB169YRFhbG8eNdLsLvEdVD74HkvLYl//oFdNDSLm9uTWPd8XzunNr3PJqiKObpridtKYcPH2bTpk3s3buXGTNmcOuttxIWFtbnclUPvQeSv5/h0vc56O0ND/UmPsiTr47m6lquoij2R0rJQw89xGuvvUZUVBS/+MUveOKJJ3QpWwX0HkgvrsHf04UAL31OF2mjpV3C2ZdRSmGlWmSkKI7s73//O1FRUcybNw+AZcuWcerUKbZv397nslVA74G8ijrC/dwsUvY148KQEtaqOemK4tCWLFnCJ5988v33RqORgwcPMmvWrD6XrQJ6D+SV1xPm626RshOCvRke6s0XR1TaRVGU3lEBvQdyK+oI97VMDx3g5sTBHDpXzv6M7ja3VBRFuZgK6Gaqbmimqr6ZUAv10AFumxxFoJcrf9lyxmJ1KIqiDUzau960UQV0M+VX1AFYLIcO4O5iZOmsOHallnAgU/XSFcUS3NzcKCkpseugLqWkpKQEN7eexRs1D91MueXa7BNL5dDb3DElmuXb0/jL5hQ+uH+KRetSlIEoMjKS7Oxs7P1MBjc3NyIjI3t0jwroZspr7aGHWTCHDlov/cGZ8by4NpmkzFISY/wtWp+iDDTOzs7ExsbauhkWoVIuZsotr0cICPGxbEAHuGNqFIFeLvxlS4rF61IUxXGogG6m/Ip6Ar1cddtlsSseLk48ODOeb1OKOZpdbvH6FEVxDCqgm8nSUxYvdFOiljv7NqXYanUqitK/mXPAxUohRKEQosMtwYQQdwghjrZ+7BZCjNO/mbaXV2G5RUUd8fNwISHYiyQ120VRFDOZ00P/J7Cwi/czgFlSyrHAC8AKHdplV6SU5JXXEWbBKYsdSYwexMGzZZhM9ju9SlEU+9FtQJdS7gA67SZKKXdLKctav90L9GyeTT9Q1dBMTWOLxWe4XOiS6EFU1jeTVlRt1XoVRemf9M6h3wes6+xNIcQSIUSSECLJ3ueAtpdnpTnoF2qbsph0tqybKxVFUXQM6EKIOWgBvdMd46WUK6SUiVLKxKCgIL2qtrhcK6wS7UhMgAcBni4kZaqArihK93RZWCSEGAu8AyySUpboUaY9sVUPXQjBxOhBHDyrBkYVRelen3voQogoYBVwl5TSIXeVyq+owyAg2Fvfgy3MkRg9iMySWoqq1HmjiqJ0zZxpi/8G9gDDhBDZQoj7hBBLhRBLWy/5XyAAeEsIcVgIkWTB9tpEbkU9wd5uOBmtP20/MWYQAIfOqbSLoihd6zblIqW8rZv37wfu161FdiivwvpTFtuMjvDFxWjg4NkyFowKtUkbFEXpH9RKUTPkldcTbuX8eRtXJyNjIn3VAiNFUbqlAno3pJTkVdQTauU56O0lRg/ieE4l9U0tNmuDoij2TwX0blTUNVHXZP1FRe1dEj2IxhYTx3IqbNYGRVHsnwro3Wg72CLczzYpF9ACOqDmoyuK0iUV0LthrYMtuhLg5UpcoKea6aIoSpdUQO9GboVtFhVdaES4D2cKqmzaBkVR7JsK6N3Ir6jDySAIssGiovYSgrw4V1qrBkYVRemUCujdyCuvJ8THDaNB2LQdQ0K8kBK186KiKJ1SAb0buRV1Ns2ftxkS7A1AaqEK6IqidEwF9G7kVdQTZsMZLm1iAj0wGoQK6IqidEoF9C5IKcmvqCfUx7b5c9BWjEb7e5BSoAK6oigd02X7XEdV3dBMQ7PJ5gOibRKCvUgpVDNd9NLY0si6jHVkVGSQXZ1NbnUu44LG8eiER/Fw9rB18xSlx1RA70JpTSMA/p72EdCHhHjxzalCGptNuDipX676oqmliZ9v+znbsrfhZHAi3DOcQPdAPkz+kG1Z23jh0hdIDE20dTMVpUdUQO9CcbUW0AO8XGzcEk1CsBfNJsnZkhqGhHjbujn9VlNLEz/frgXzZ6Y8w81Db8ZoMAKQlJ/Ec7ue454N93D3yLt5IvEJhLDtDCdFMZfq5nWhrYce4GkfAb1tpkuKGhjttSZTE7/Y8Qu2Zm3lmSnPcNvw274P5gCJoYl8vvhzbhx6I++dfI/Vqatt11hF6SFzDrhYKYQoFEIc7+R9IYT4qxAiVQhxVAgxUf9m2kZJtXZKUICXfaRc4oO8EEJNXeyLF/e+yJZzW3hq8lPcNrzjrf49nD14bupzTAyeyKtJr1JcV2zlVipK75jTQ/8nsLCL9xcBQ1o/lgBv971Z9qHEznro7i5GIvzcVQ+9lw7kH+DzlM+5Z/Q93DHiji6vNQgDz09/nrrmOl7a/5KVWqgofdNtQJdS7gC6Ol3hWuA9qdkL+AkhwvRqoC2V1jTi6WLEzdnY/cVWMiTYixRH3tOlsRbqynUvtsnUxIt7XyTcM5yHxj1k1j1xvnE8OPZBNmRuYFvWNt3bpCh602NQNALIavd9dutreTqUbVMl1Q3428mAaJshId7sSiuhxSRtvh2BLhqq4ND7kH0A8o9BaRpIE/hEQsgoCBsLifeCT3ifqvng5AekVaTx+tzXcXcyf6HYvaPvZX3men6393ckhiTi5eLVp3YoiiXpMSjaUVSRHV4oxBIhRJIQIqmoqEiHqi2rpKbRbqYstkkI9qKx2URWaa2tm9I3zQ2wdzn8ZTxseBpykiBwKMz8BVzxa4ieDhVZ8O2f4PVE2PEqNNX3qqr8mnzePvI2swfPZvbg2T2619nozK+n/5rC2kJWHl/Zq/oVxVr06KFnA4PbfR8J5HZ0oZRyBbACIDExscOgb09KqhvtYh+X9hKCtR5iSmE1MYGeNm5NL2V8C2uWQfk5iLkMrvgNRF7S8bWlGbDxV/DNC/Dd+3D1axA/p0fVvbz/ZaSUPDX5qV41d1zQOOYMnsNnZz5jydgluDnZ178JRWmjRw/9C+Du1tkuU4EKKWW/T7eAlkP3t5MB0TY/BPR+mkdP/go+uAGc3ODOVfDjLzsP5gD+sXDrh3D3GjC6woc3won/ml3d4cLDbD63mQfGPkCEV0Svm33nyDspbyjnq/Svel2GoliaOdMW/w3sAYYJIbKFEPcJIZYKIZa2XrIWSAdSgb8DyyzWWiuSUlJS02A3Uxbb+Lg5E+rjRmp/3NPluw/h07sgbBzcuwESLgdzF+3EzYb7N0HkJPjPvXD4I7Nue/f4u/i6+nLniDt7324gMSSR4f7D+eDkB0hp979cKgNUtykXKWXHk3V/eF8CD+vWIjtR1dBMU4u0mymL7SUEe5Ha3/ZF3/c3WPdLiJsDt3wArr0YXHTzhTs/h49vh9UPQVMtTLq/08vTK9LZmrWVJWOX9HlvFiEEd428i2d3PsuevD1MD5/ep/IUxRLUStFOlFa37eNipwG9sBqTqZ/0FDN2wLonYfjVcPsnvQvmbVw84bZPYOgi+PrnkPxlp5e+d+I9XIwunS4g6qmFMQsJcAvg/ZPv61KeouhNBfROlNS0rRK1v4A+JMSL2sYWclsPsLZrtaWw6kEIiIcbVoCTDiksZze4+V8QPhFWL4OStIsuKaot4ou0L7gu4ToC3AP6XifgYnThluG3sDNnJ+kV6bqUqSh6UgG9EyVtG3PZ2bRFgLhArYebXlRj45Z0Q0r44lGoKYIf/UPrXevFyVUL6sIAn/4Yms7/4fZh8oe0yBbuHnm3fnUCNw+9GReDCx8lm5fDVxRrUgG9E98v+7fDHnp8kBYYM4rtPKAf+hec+gqueB7Cx+tfvl8U3PB3KDgGa5/4/uXqxmo+Pf0pV0RdQZRPlK5VBrgHsCh2EV+mfUl9c+/mxSuKpaiA3okf9kK3v4Ae5O2Kp4uRdHseGC1Jg3VPaYOgUy04Zj50vrYY6bsP4PC/AViVsoqqpiruGX2PRaq8Mu5Kaptr2Z272yLlK0pvqYDeieLqBrxcnexqH5c2QgjigrxIt+ce+pbfgMEI1y8Hg4X/mc1+GqKmwYankdXFfHbmM8YFjWN04GiLVDcpdBK+rr5sOrvJIuUrSm+pgN4Je1xU1F5soKf9plxyv4OTa2DaI+Adavn6DEa4+s/QUMWRDf9DZmUmPxryI4tV52xwZu7guWzL2kZjS6PF6lGUnlIBvRMl1fYd0OOCPMkpr6O+qcXWTbnYlt+Cuz9Ms+LyhOARMO1hVuXtwt3oyvyY+Ratbl70PKqbqtmbt9ei9ShKT6iA3omSmkYC7XBAtE1soCdSwtkSO9ukK+NbSPsGLvsZuPlYteqa6Y+w3suThfUteBos+3c3NWwq3s7ebMjcYNF6FKUnVEDvRGlNg1330OOD2qYu2tHAqJRa7tw7vMsVnJayIXcndUJwQ1E27H3LonU5G52ZEzWHrVlbaWppsmhdimIuFdA7IKWktKbR7vZxaa9tp0W7Ghg9vU7b13z2k+Bs/p7jelmVsopY31jGRc2FbS9DtWW3aJ4fPZ+qxir25e+zaD2KYi4V0DtQWW+/+7i08XJ1Itjb1X4GRqWE7S+DfzyM79tGWL2RXp7OkaIj3JBwA2LB76C5Dnb+yaJ1Tgufhqezp5rtotgNFdA78MPh0PYb0EEbGLWblEt2EuQdhqkPgVGPbfZ7ZlXKKpyEE9fEXwOBQ2Dc7XDgH1CRY7E6XYwuzB48my3nttBkUmkXxfZUQO/AD4uK7DflAhAb6GU/PfQD74CLN4y71epVN5ma+DL9S2ZGzvxh35ZZv9SOstvxR4vWPS96HhUNFRwsOGjRehTFHCqgd6D4+31c7LuHHh/kSVltE2U1Np4LXVMMJ1ZpwdzV2+rV78vbR2l9KYsTFv/w4qBouOQn2grSUsttpDUtbBrOBmd256hVo4rtqYDegVI73selvVh7GRg99B60NNpkZgvA2vS1eDt7c1nEZee/MfMJMDjDtpcsVreHswfjg8erbQAUu2BWQBdCLBRCnBZCpAohLjqYUQjhK4T4UghxRAhxQghhmU00rKS0detce562CBBnD1MXTS2Q9K52NmjwcKtXX9dcx5ZzW5gXMw8X4wV/X96hMPkBOPopFCZbrA3TwqZxuuw0xXXFFqtDUcxhzhF0RuBNYBEwErhNCDHygsseBk5KKccBs4H/E0LYdzTsQnF1I96uTrg62d8+Lu1FDnLHySBsm0c/swEqzmmB0wZ2ZO+gtrmWK2Ov7PiCGY9r2/Z+a7kZL9PCpwFa6kdRbMmcHvpkIFVKmS6lbAQ+Bq694BoJeAshBOAFlALNurbUikprGvG383QLgLPRQFSAh20D+oG/awuJhl1lk+rXpq8lyD2IxJDEji/w8Ndy6cc/h7JMi7RhhP8IfFx82JO7xyLlK4q5zAnoEUBWu++zW19r7w1gBJALHAMek1KaLixICLFECJEkhEgqKrLsoo++KLHzVaLtxQV62u6gi9IMbZl/4j02mapY0VDBtznfsiBmAUZDF79NTXtE28Br118t0g6jwciUsCnsydujDpBWbMqcgN7RsewX/qtdABwGwoHxwBtCiIs28pBSrpBSJkopE4OCgnrYVOspqW60y5OKOhIX5EVGSY1tzhc9uUb7bIOpisD387+viuvmtwOfMBh3mzbjparAIm2ZFj6NwtpCMioyLFK+opjDnICeDQxu930kWk+8vXuAVVKTCmQA1h8h00lJTaPdT1lsExvoSWOziZxyG5wvenK1dq6nn76nAplrbfpaoryjGBUwqvuLL30MTE0W2+NlWpiWR9+Tp9Iuiu2YE9APAEOEELGtA523Al9ccM054HIAIUQIMAzol6fomkySsppGu5+y2CYu0EbH0ZWd1fY9H3nhcIp1FNYWsj9/P1fGXYk2dNONgHgYeZ22erSuXPf2RHpHEuUdpfLoik11G9CllM3AI8AGIBn4VEp5QgixVAixtPWyF4DpQohjwBbgSSllv5zDVVnfRLNJ9pscemzr+aJWn7qY3Poz3UYBfWPmRiSSRbGLzL9pxuPQWKWtarWAaeHTOJB/QG0DoNiMWfPQpZRrpZRDpZTxUsoXW19bLqVc3vp1rpRyvpRyjJRytJTyA0s22pLaDocOtOOdFtsL8nLF29XJ+ouLTq6BsHHgH2vdelutz1zP0EFDifONM/+msLGQMA/2LYcm/Q94nhY2jdrmWo4WHdW9bEUxh1opegF7Phy6I0II4oO9SCmwYg+9IlvbJtdGvfO86jyOFB1hYczCnt88/VGoKYJjn+rerklhkzAIg1o1qtiMCugXaNtpsb8EdIChIV6kFFoxoJ9sTbeMsFG65exGgN4dMxc7E0LHwO43wHTRzNo+8XHxYXTgaHUsnWIzKqBfoG1jrv6ScgEYEuxNcXWD9TbpOrkGQkZDYIJ16rvAxsyNjPAfQbRPdM9vFgKmPQrFpyF1s+5tmxI6hRPFJ6hutJNtjZUBRQX0CxRW1mMQ2PV5ohdKCNH2dEm1xsBoZS5k7dVmjNhATnUOR4uP9u0Q6NE3aKtb97yuX8NaTQ2bSotsIakgSfeyFaU7KqBfIL+ynkAvV5yM/efRDAnWAvqZgirLV5b8lfbZRvnzTZna6UALYhb0vhCjM0xdChk7IE/fAcxxweNwNbqqfV0Um+g/UctKCiobCPFxs3UzeiTCzx1PF6N1BkbPrIeAIRA01PJ1dWBD5gZGBYxisPfg7i/uysQfg4sX7HlDn4a1cjW6MiF4gsqjKzahAvoFCirrCfHpP/lz0Ga6JAR7kWrpgdGmeji7CxKusGw9nciuyuZ4yfG+9c7buPvBxLu1Tbt0PqZuStgUUstT1Xa6itWpgH4BLaD3rx46QEKwNymFFk65nNsNzfUQP9ey9XRiQ+YGoJezWzoyZal2TN2Bv+tTXqupYVMB2J+3X9dyFaU7KqC309DcQlltU78M6ENCvCiobKCizoKrFNO+0U4AirnUcnV0YUPmBsYEjiHC68LNPntpUDQMv1o7oKOxVp8y0bbT9Xb2Zl++yqMr1qUCejuFldoc9ND+GNBbB0ZTLdlLT9sKUVO1AyOsLKsyi+TSZH3SLe1NXQb15XD0Y92KNBqMJIYmqoFRxepUQG+noFJbDh7cz3LoAENDtMOZLTYwWpUPBcch4XLLlN+NDWe1dMu86Hn6Fhw1FcLGw963dV1oNCVsCjnVOWRVZXV/saLoRAX0dvJbA3qob//roUf4uePmbLDcitG0rdpnG+XPN2ZuZGzgWMK9wvUtWAitl158BtK/0a1YlUdXbEEF9HYKWlMuId79L6AbDNpMF4vNRU/7BjwCIWSMZcrvwrnKcySXJus3GHqhUdeDVyjs0W+v9DjfOILcg1TaRbEqFdDbKaysx8XJgJ+Hs62b0itDgr0tM3XRZIL0rRA/BwzW/yfz/d4t0RYK6E4uMPl+SNsChad0KVIIweSwyezL36eOpVOsRgX0dvJb56CbdWCCHRoS4kVeRT1V9TrPdCk4ru1QGG+b/PnGzI2MDRpLmFeY5Sq55B5wctO21tXJlNAplNaXcqbsjG5lKkpXzAroQoiFQojTQohUIcRTnVwzWwhxWAhxQgixXd9mWkdBZX2/TLe0GRKsDYzq3ktP26J9jp+jb7lm+D7dYqneeRvPQBhzExz5GGpLdSlyevh0AHbl7tKlPEXpTrcBXQhhBN4EFgEjgduEECMvuMYPeAtYLKUcBdykf1Mtr6CygZB+OCDapm3qou4Do2nfQPAo8A7Vt1wzWDzd0t7Uh6C5Dg69p0txIZ4hDBk0hF05KqAr1mFOD30ykCqlTJdSNgIfAxfuzHQ72iHR5wCklIX6NtPypJT9voc+2N8DFycDKXoOjDbWwrm9Numdg7aYaFzQOMumW9qEjIKYy2D/36GlWZciZ4TP4FDhIWqarHyilDIgmRPQI4D2k2mzW19rbygwSAixTQhxUAhxt14NtJbqhmZqG1sI9e1/c9DbGA2C+CCdD7s4txtaGm0S0M9WnuVU6Snr9M7bTH0IKrPh1Je6FDcjYgbNpmY1fVGxCnMCekcjhBcO2zsBlwBXAQuA54QQF23HJ4RYIoRIEkIkFRUV9bixltS2qKg/Lvtvb2iIzsfRpW8DowtETdevTDOtz1gP6Lh3izmGLgS/aNirz+DohOAJuDu5qzy6YhXmBPRsoP1epZFAbgfXrJdS1kgpi4EdwLgLC5JSrpBSJkopE4OCgnrbZotom4Me3I9TLqCtGM0pr9NvT5f0bTB4Crh46FNeD6zPXM+E4AmEeloxd28wwpQHtUM8cr/rc3HORmemhE1hZ85ONX1RsThzAvoBYIgQIlYI4QLcCnxxwTVrgMuEEE5CCA9gCpCsb1MtK7+i/64SbW90hC8AJ3Iq+l5YTTHkH4O42X0vq4fSytNILU/Vf+8Wc0y4U9srXade+ozwGeRU53C28qwu5SlKZ7oN6FLKZuARYANakP5USnlCCLFUCLG09ZpkYD1wFNgPvCOlPG65ZuuvoKot5dJ/c+gAY1oD+jE9AnpG6+zTOOvnzzdkbkAgrJs/b+PmC+Nv1/ZKryroc3HTI9T0RcU6zJqHLqVcK6UcKqWMl1K+2Pracinl8nbXvCKlHCmlHC2lfM1C7bWYgop6vN2c8HBxsnVT+sTf04UIP3d9Anr6NnD1hfDxfS+rB6SUrM9cT2JoIkEeNkrNTVkKpmY48E6fixrsPZgYnxh25uzUoWGK0jm1UrRVfzx6rjNjInz7HtClhLRtEHuZlle2ojNlZ8ioyGBhzEKr1nuegHhtgDTpH9pJTX10acSlJOUnUd/c97IUpTMqoLcqqKrvl/ugd2RMpC9nS2r7NjBamg4V52ySP9+QuQGDMHBFtG2OuvvetGVQWwLHPu1zUZeGX0p9Sz2HCg7p0DBF6ZgK6K0KKur75T7oHdFlYDR9m/bZyvnztnTL5NDJ+Lv5W7Xui8RcBiGjtV0Y+zhDJTE0EReDC9/mfKtT4xTlYiqgAyaTpLDKsVIu0MeB0fRt4DtYSz1Y0cnSk2RVZdk23dKmba/0ouQffsD1kruTO5PDJrM9e7uavqhYjAroQElNI80m6TApl7aB0aO9DeimFsjYAXGztKBmRRsyNuAknGyfbmkz5kbwDIa9fd8rfXbkbLKqssioyNChYYpyMRXQab9K1DFSLqD10o/3NqDnHdHO2bRyusUkTazLXMe08Gn4uvpate5OObnCpPshZSMU9W0b3FmDZwGwLXubDg1TlIupgI7jLPtvr08Do2mtR7HFztS3Ud04WHCQ/Jp8roq7yqr1divxXjC69rmXHuoZynD/4WzP6pe7Syv9gArotDt6zpECel8GRlM3Q9g48ArWuVVd+zr9a9yd3Jkz2DY7O3bKKwjG3QqHP4Lqvu1BNCtyFoeLDlNeX65P2xSlHRXQ0XroQkCQt2OlXICe59HryiBrHwyx7grNxpZGNp7dyNyouXg4W3/fmG5Nf1TbdXL/ij4VM3vwbEzSpGa7KBahAjpaQA/wdMXZ6DiPY1BvV4ymbQVpgoR5lmlYJ3bm7KSqsYqrYu0s3dImcAgMuxIO/B0ae7+3+ciAkQS6B7Ita5tuTVOUNo4TwfqgoPUsUUfTq4HR1M3g5geRiRZpU2e+Tv8afzd/poVPs2q9PXLpT7XfYL77sNdFGISBWZGz2JW7i6YWnc9+VQY8FdCBvArHWSXa3vcDo7VmBg6TCVI2QcLlVl3uX91Yzfbs7SyIWYCTwY730omaCpGTYc8bfTrRaFbkLGqaajhQcEDHximKCui0mCQZxTXEBXnauim6a8ujH881s5eefxRqCq2ebtl8bjMNLQ32N7ulI9MfhfKzkHzhDtLmmxo+FVejq5rtouhuwAf07LJaGppNDAn2tnVTdDc2UgvoB8+WmXdD6ibtc8LlFmpRx75O/5rB3oMZGzjWqvX2yvCrwD8Odv+119sBuDu5MyVsilo1quhuwAf0tuPaEkK8bNwS/fl5uDAq3IddqcXm3ZCyGcInWHW6YmFtIfvz93Nl7JUIK69K7RWDUeul534Hmb3fDnfO4DnkVOdwuuy0jo1TBjoV0FsPVE4IdryADjAjIZBD58qobewm51tbCtn7rZ5u+TLtS0zSxOL4xVatt0/G3QYegbDrL70uYm7UXAzCwKazm3RsmDLQmRXQhRALhRCnhRCpQoinurhukhCiRQhxo35NtKyUwipCfdzwcXO2dVMsYnpCIE0tkgOZ3aRd0lunKw6xXkCXUrI6dTUTgycS5RNltXr7zNldOwAjdRMUnOhVEf5u/lwScgmbz27WuXHKQNZtQBdCGIE3gUXASOA2IcTITq57Ge2oun4jtbCaIQ6YbmkzKWYQLkZD92mXlM3gPggiLrFOw4AjRUfIrMzkuoTrrFanbibdB84esPv1XhcxL3oe6RXppJWn6dgwZSAzp4c+GUiVUqZLKRuBj4FrO7juUeBzoFDH9lmUySRJLax22HQLgIeLExOj/diZ0kVAb2mGlA2QcIVVpyuuSVuDu5M782NscG5oX3n4w8Qfw7HPoCK7V0VcHqUNPqu0i6IXcwJ6BJDV7vvs1te+J4SIAK4HujwmXQixRAiRJIRIKirq254YesitqKO2scUhZ7i0NyMhkJN5lZRUN3R8QcZ27WSekddZrU11zXWsz1jPvOh5eDr30ymj05ZpM132vt2r24M9ghkfNF6lXRTdmBPQO5p6cOFcq9eAJ6WULV0VJKVcIaVMlFImBgXZ6PDfdtoGRB055QJwaUIgAHvSSzq+4MQqcPXReuhWsuXcFqqbqvtnuqWNXxSM/hEc/Ke2grQXroi+gtNlp8mqzOr+YkXphjkBPRsY3O77SCD3gmsSgY+FEJnAjcBbQojr9GigJaW2TVkMcuyAPibCF283p47z6M2NkPyltk+Js/VWy65JXUOEVwSXhFgvZ28Rl/4UGqvhwD96dXvbQR6bzqm0i9J35gT0A8AQIUSsEMIFuBU4b5mclDJWShkjpYwB/gMsk1Ku1ruxeksprCLQy5VBni62bopFORkNTI0LYGdHAT19K9RXwKjrrdaevOo89uXt49r4azGIfj5zNnQMxF8O+/4GTfU9vj3CK4JRAaNU2kXRRbf/N0kpm4FH0GavJAOfSilPCCGWCiGWWrqBlpRSWM0QBx4QbW9GQiBZpXWcK6k9/43jq8DNF+LnWq0tq9NWI5FcE3+N1eq0qEsf07ZMOPLvXt1+RfQVHCs+Rl51ns4NUwYas7pHUsq1UsqhUsp4KeWLra8tl1JeNAgqpfyJlPI/ejdUb1JKUgsce8pie2159PN66U31cHotDL8GnKzzW0qLqYVVKauYGjaVSO9Iq9RpcbEzIWy8NoXR1OUwUofmRWtz/zee3ahzw5SBpp//vtt7BZUNVDU0D5geenyQJ6E+bufn0dO2QEMljLZeumVX7i7ya/K5aehNVqvT4oSAGf8DpWlw6qse3x7tE83ogNF8mfal/m1TBhQ73qvUslIKqwBI0GnKYlFtEYeLDlNcV0xxXTGNLY3cPOxmBnsP7v5mKxBCMGtoEF8fy6OusQV3F6OWbnH3h9hZVmvHZ6c/I9A9kDlRdnbMXF+NWAyDYmHna9rXPdyX5pr4a/jD/j9wuvQ0w/yHWaaNisMbsD30tk25+ppykVLyZdqXXLP6Gn627Wf8ft/veefYO3xw8gOuW30dbx5+k7rmOj2a3GfXTginuqGZTckF0FgLp9fByMVgtM62B/k1+ezI2cH1CdfjbHCwrRa+37TrEJzd1ePbF8UuwsngxBdpvd+WV1EGbkAvrGaQhzMBfZjhUtFQwS93/JJndj7DsEHD+PDKD9l681YO3XmI9T9az+XRl7P8yHKuW30d+/L26dj63pkaG0C4rxv/PZStTVVsqoFRN1it/lUpq5BS8qOhP7JanVY1/nZt065v/9TjWwe5DWJmxEy+Tv+aZlPvD89QBrYBG9BTC6sYEuzd6y1bz1We46Yvb2Lz2c38dMJPWblgJWODxhLoHojRYCTEM4Q/zvwjKxesxNXJlWWbl7Ezp/fbrerBYBBcNyGCHSlFNO1+CwKHagN6VtBsaubzlM+ZHjGdCK+I7m/oj5zdYfoj2thE1v4e3744fjEl9SXsyd1jgcYpA8GADOhSSs4UVPd6D/T8mnwe2PgA9c31vLfoPR4Y+wDGTvZAmRQ6ifcXvU+8Xzw//eanfJtt29Peb5gYwViZgnPBYZi8pMe53t76NvtbCmsLHWswtCOTHgCPANj2hx7felnkZfi6+qrBUaXXBmRAL65upKKuqVczXErrS1myaQmVjZUsn7ecMUFjur3H19WXv8//Owl+CTy29TGbBvWEYG/+x3sL1cJT29fbSj478xnB7sHMirTeAKxNuHpp89LTvoFzPUuzuRhdWBSziG+yvqGqscpCDVQc2YAM6Ml5lQAMDenZDJeqxiqWblpKbnUur899nZEBF+0i3KkLg/runN09qls3lblc1rSbfzfN4ky5dY4/y6rMYmfOTm4YeoN9HwKtl0n3g2cQbPt9j29dHL+YhpYGNmaqOelKzw3IgH4gsxSjQTBusF+P7ntp/0uklKXwp9l/IjE0scf1tgX1ON84Htv6GEn5ST0uo8+SViJkCx+a5rPqUI5Vqvzo1EcYhZGbh95slfpszsVT66Wnb4OzPcuHjw4cTYxPDKtTV1ukaYpjG5ABfV9GKaPCffByNb+3eLr0NF+mfcndo+5mZmTvBxJ9XX3527y/EeYVxsNbHuZo0dFel9VjTfWQ9C5i2CLiho5mzeEcTCbL9tJrmmpYnbqa+THzCfKw/Q6bVpN4H3gG97iXLoTgxqE3crjoMCdLTlqocYqjGnABvaG5hcNZ5UyO8e/RfX8++Gd8XH24b8x9fW5DgHsA78x/hwD3AJZuXmq9/3FPrILaYpiylOsnRJBXUc/2M5bdl/6LtC+obqrmjhF3WLQeu+PiATMeh4wdkLqlR7deP+R63J3c+eDkBxZqnOKoBlxAP5pdQWOzicmx5gf0vXl72ZW7iwfGPICPi48u7Qj2COad+e/g5ezF/Rvvt3xPvaUJvv0/CB4JsTNZMCqUcF833tqWarEqTdLER8kfMSZwDGODxlqsHrs16T5t9ej6p7XnbyYfFx+uT7iedZnrKKq1/UEwSv8x4AL6/oxSACaZ2UM3SRN/PvhnwjzDuHX4rbq2JdwrnH8u/Cd+rn48sPEBy+bUk1ZCSSpc/jwIgYuTgSUz4ziQWfb9M9Hbntw9ZFZmcvuI2y1Svt1zcoUFv4fi03DgnR7deseIO2gxtfDJ6U8s1DjFEQ24gL4vo5ShIV5m74G+IXMDJ0tO8uiER3E1uurennCvcN5d8C4hniE8tPkhdudaYPZLXTlse0lbRDR0wfcv3zIpigBPF97Yaple+gfJHxDoHsiC6AXdX+yohi3Stibe+geo6eag7naifKKYNXgWn57+lIaWTo4OVJQLDKiA3txi4tDZMrPTLS2mFl7/7nWGDRrGVXFXWaxdIZ4hvLvgXaJ8onh488N8evpTfSv49lXtiLT5L563kMjdxch9l8Wy40wRx7IrdK0ysyKTnTk7uXnozThbaa8YuyQELHxJ22Zhy297dOtdI+6irKGMtelrLdQ4xdGYFdCFEAuFEKeFEKlCiKc6eP8OIcTR1o/dQohx+je175LzqqhuaDY73bI7dzdZVVk8MPYBi5+sE+AewLsL32Vq+FRe2PsCv9v7O5p6kHftVFmmdprO+Nsh7OI89l1To/F2c+JNnXvp7518D2eDMzcNc/CVoeYIGgaTH4RD70HuYbNvmxQ6iaGDhvJ+8vtIaZ01A0r/1m2UEkIYgTeBRcBI4DYhxIUrajKAWVLKscALwAq9G6qHfRnaIcnm9tBXpazC382fuYOtc5qPj4sPb8x9g3tH38snpz/hgU0PUNlY2bdCN/8aDE4w91cdvu3t5sxPpsew/kQ+KQX6rE4sritmTeoaFscvJtA9UJcy+73ZT4JnIHz1P9Bi3uZbQgjuGnkXKWUpNt8HSOkfzOl2TgZSpZTpUspG4GPg2vYXSCl3Synbjj3fi3aQtN05kFlKlL8HYb7u3V5bXFfMtqxtLI5fbNWUgdFg5PFLHufly17mSNERntzxJC29OAUHgNPr4cR/tW1dfcI7veyeS2Nxdzbqlkv/96l/02Rq4sejfqxLeQ7BzReufAVyv4PdfzH7tqtiryLCK4LXv3sdkzRZsIGKIzAnoEcAWe2+z259rTP3Aes6ekMIsUQIkSSESCoqsu50LCkl+zNKzU63fJH2Bc2ymeuHWO80n/aujLuSpyc/zc6cnbx15K2eF1CZC6sfgtCxcNnPu7zU39OFn1waw5rDuRzP6Vsuvbaplo9PfczcqLnE+sb2qSyHM+p6GHmdNkBdYN7aA2ejM8vGLyO5NFkdJK10y5yA3tF2fB0m9IQQc9AC+pMdvS+lXCGlTJRSJgYFWXfVYGphNWW1TUwxI90ipWRVyiomBk8kzjfOCq3r2E1Db+KGITew4ugKtpzrweIUUwt8/gA0N8CN72rT57qxdFY8fh7OvLz+VB9arKWpKhsruWf0PX0qx2Fd+Sq4esOaZWanXq6KvYp433he/+51tVe60iVzAno20P4ctUgg98KLhBBjgXeAa6WUJfo0Tz/7M7W51ubkz5MKkjhbedbmBzEIIXhmyjOMDhjNszufJb0i3bwbd7wKZ3fCVa9CYIJZt/i6O/PInAS+TSnm25Te/fbUZGrivZPvMTF4IuOC7HJc3Pa8guCq/+tR6sVoMPLohEfJrMxUW+sqXTInoB8AhgghYoUQLsCtwHnnZAkhooBVwF1SyjP6N7PvdqeWEOztSnSAR7fXfp7yOd7O3t+fxm5LrkZX/jznz7gaXXli+xPdz3zJ+Ba2vwRjb+nx9rh3TYsmcpA7L6071as9XjZmbiSvJk/1zrsz6noYea02Nz3noFm3zI2ay+iA0bx95G0aWxot3EClv+o2oEspm4FHgA1AMvCplPKEEGKpEGJp62X/CwQAbwkhDgshbLCNYOeqG5rZcqqA+aNCuj2hqKKhgk2Zm7gy7krcnbofPLWGUM9QXrj0BVLKUlh+dHnnFxachE/uAP94rRfYw8MrXJ2MPDF/GCdyK/niyEW/hHXJJE384/g/iPON69PmZQPG1a+Bdyh89hNtjUA3hBD8dOJPyavJU6tHlU6ZNblaSrlWSjlUShkvpXyx9bXlUsrlrV/fL6UcJKUc3/rR871lLWjjiXzqm0xcN777o89Wp66m0dTIj4bY17mXMyNnsjh+Mf849g9OlJy4+ILyc/DBDeDsAXet0vK0vbB4XDijwn14ZcNp6pvMn12z5dwWUspSWDJ2icXn7DsED3+46Z9QmQerl4EZ88ynhU9jWtg03jr8FgU1BZZvo9LvDIj/81YfziVykDuXRA/q8rqGlgb+deJfTA6dzIiAEVZqnfl+OemXBLgF8Kudvzr/1+6aEnj/BmishTs/B7+oXtdhMAievXIEOeV1vLUtzax7TNLE8iPLifGJYWHMwl7XPeBEJsL8F+D0Wtjzhlm3/Grqr2gyNfHivhfVYiPlIg4f0IuqGtiVWsy148O7TbesSV1DUV0RD4x9wEqt6xlfV1+en/48qeWpLD/SmnqpK4OPbtJ66Ld/DCGj+lzP9IRArh0fzvJtaaQVVXd7/TfnvuFM2RkeHPdgp2erKp2YshRGLIZNz0Pmrm4vj/KJ4uHxD7M1aysbz6pTjZTzOXxA//poLi0mybXdpFuaTE2sPL6SsYFjmRI6xUqt67m21MvK4ys5fnY7/PNqyD8GN/8LoqfrVs+zV43A1dnAc6uPd9kTNEkTbx95mxifGBbFLNKt/gFDCLj2DQiIh49vg8Lkbm+5a+RdjAwYye/3/Z6KBn334FH6N4cP6KsP5zIizKfb80PXZawjpzqHB8Y+0G1P3taenPwkga6DeGbLI9SXZsDtn2q7+uko2NuNXy4czu60EtYc7nyAdOu5rZwpO8OSsUtU77y33Hy1VJmTO3zwI6jI7vJyJ4MTv5n+GyoaKnjlwCtWaqTSHzh0QM8sruFwVjnXje982Ttovcx3jr3DkEFD+sUMDZ/Sc/y2oIAMI/x12q0QP8ci9dwxOYrxg/343dcnqai9eLrkeb3zWNU77xO/KLjzP9BQpQX1bma+DPcfzj2j72FN2hp1/qjyPYcO6F8cyUUIWNxNQN98djMZFRk8MMbyuyr22ZFP4J0rmN7Ywi2D5/FB1mYO5B+wSFUGg+DF60dTWtPILz8/ctHc9K/Tv+Z02WmWjF2Ck8H881mVToSOgVs/hNJ0+PAmqO364JFl45cxJWwKv9n9G/bl7bNSIxV7ZufRq/eklKw+nMPkGP8uN+Oqa67jzcNvEu0Tzfzo+VZsYQ81N8LaX8B/l0DERHhwBz+77HdEekfy3K7nqGmqsUi1o8J9efaqkWw4UcBL7bYFqG6s5v+S/o8xgWMsulf8gBM7U9uuIe8IvLuoy/SLs8GZP83+E9E+0Ty+9XHSy81cSaw4LIcN6F8fyyO9qIYbL+l648cX975IRkUGz0x+xn5zwAUnYOV82L8Cpj0Cd68B7xA8nD14ccaL5Fbn8sLeFyw2je3eS2O4e1o0K3ak8+G+swC8feRtSutLeXbKs/b/W01/M+JquHOVtsHaP+ZD0elOL/Vx8eHNK97ExejCsi3LKK4z/1QkxfE45P+JVfVN/PbLk4wK9+GGiZ0H9P+m/Jc1aWtYMnYJ0yP0myGim+YG+OZ38LeZ2rTEm9+DBS9Cu+18JwRPYNn4ZXyd/jUfJn9okWYIIfjfq0cyZ1gQ/7vmBB8f3s+HyR/yo6E/YlRg36dJKh2IvQzuWQumZli5AFI732kxwiuCNy5/g5K6Eu5ceyenSzv/AaA4NocM6K9tTqGouoHfXTcao6HjGSunS0/z4r4XmRI2hYfGPWTlFnajpRmO/QeWz4Adr8DoG+HhA9r+Hx1YMnYJcwfP5dWkV9mft98iTXIyGnj99okMCfHihT2/x1m48+j4Ry1Sl9IqdAzcuwG8w7SB0k3PQyd7+YwOHM3KBStpamnirnV3sensJis3VrEHDhfQT+ZW8s/dmdw2OYoJUR2vDK1oqOCJ7U/g4+LDS5e9ZD+ploYq2P93eH0ifH6fthz8js/hhr+BZ0CntxmEgRdnvEiUTxRPbH+C3Oqe7cNiLi9XJ+6dX4bBI43ynCt46rN0ymvVRlEW5R8L92+BS34Cu16Dd6+EsrMdXjomaAwfX/0xQwYN4WfbfsZfDv1FbeQ1wAhbLR9OTEyUSUn67uFlMklu+tseMopr+Obns/DzcLnomvyafB7c9CBZVVmsmLeCxFAbbjvT0gSFJyF9G6RsgnN7wdQEkZPg0v+BYVeCwfyfuRkVGdz+9e1EekeyYt4KBrl1vdVBT6WVp3HH2juI9Ylljs9veWVDCoFerjx39UgWjQ61+/n7/d7xz+GLx7Q0zMyfw/SfdrjXfWNLI7/f93s+T/mcGJ8Ynpv6HJPDJtugwYolCCEOdrZflkMF9Le2pfLH9ad55cax3JQ4+KL3U8tSeXDzg9Q21fLXuX9lUugkXevvlMkEFeegOEUb4CpKhryjUHQK2npQIaMh4XIYfrUW0HsZHHfl7OKxrY8R5hnG3+b9jXCvrqdsmquioYLbv76dmqYaPr76Y0I9QzmWXcHPPzvMmYJqxkb68ssFw5kxRJ0halHlWbDhGUj+AvzjYMEfYOiCDv+97MzZyYt7XyS7Opur467msYmPEeoZaoNGK3py+IBuMkleWn+KFTvSuWpMGK/fNgHDBbnzPbl7+Pn2n+NudOetK95imP8wXeo+j5RQmaMF64ITWsAuPg3FqdBc98N1HgHa0XChYyBsnLZkv4szP3vqUMEhHvnmEdyMbrx9xdt9/rO2mFp4eMvD7Mvfx8oFK5kQPKHde5L/fpfDnzedIae8jskx/tw7I4YrRoTgZHS4jJ79SN0C656EkhStMzB1GYy58aIee31zPe8ce4eVx1cCcH3C9dw35j7dftAr1ufQAb2x2cQv/nOENYdzuXtaNM9fM+q8gdCSuhL+dPBPfJH2BbG+sSy/Yrk+/5ibG6E0TdtHJf+o9jnvKNS1WwziFwWBwyBoGAQO/eGji3y4XlLKUli6eSl1TXX8evqvmRc9r1cpESklrya9ynsn3+P5ac9z49AbO7yuobmFj/ad451vM8gpryPCz507pkZx7fgIIvzsY195h9PcCEc/gb1vaak7zyAYc5OWqouaBsYfFnvlVuey8vhKVqWsQkrJgtgFXB13NVPDpqpFYf2MQwZ0k0my7Uwhb3yTyqFz5fxiwTCWzY7/PmjVNtXyRdoXvP7d69Q213LPqHt4YOwDPT+0QkptymDBccg/DoUnoPCUFszbznc0ukLISK3HHTpW+wgZBa5evf7z6SGvOo/Htj5Gcmkyl0ZcyrOTn2Wwz8WpqM6U1pfy/O7n2Za1jVuG3cKvpv6q23taTJLNyQX8c1cme9K1kwjHD/bj6rFhzBkeTFygp8q1601KbRxm/wqt597SAG5+ED9XC+xRUyB4FBidyK/J593j7/Jl+pdUNVbh7+bPvOh5TA2byvjg8QS6q5SZvetzQBdCLAT+AhiBd6SUL13wvmh9/0qgFviJlPJQV2X2NqBX1jfxWVI27+3J5GxJLcHerjxz5QiumxBBXXMdB/IP8FX6V2zL2kZdcx2TQifxqym/Is7PjMOeW5qgLFPr7eR+98NHfduOdkKbdRA0AoKHQ9Bw7dfdwCHnzQ23J82mZj4+9TFvHH6DppYm7hhxB4vjF5MwqOuzRnfn7ObZXc9S0VDB45c8zh0j7ujxAqLM4hq+PpbH2mN5nMitBCDQy4VJMf5cEj2IoSHeJAR7EebrpoK8XhqqIX0rnFqrBfmq1hlPLl7av9XQ0RAymsaAeHY2FPF1/m52ZO+gvqUegGifaEb4jyDOL4443zhifGII8wrD29lb/R3ZiT4FdCGEETgDzEM7MPoAcJuU8mS7a64EHkUL6FOAv0gpu9yDtrcB/aP9Kfzq6x0MDZdMH+pCWEATmZXpnCw5SXpFOiZpwsfFhwUxC7gqegET/YYimmqhqRbqK6G+DOrKtc2PqguhOh+qCqAsQ9tDo63XbXDSetnhE7Q8d8gYrRfu4tnjNtuDwtpCXk16lY2ZG2mRLQwZNIQF0QuI9okmwD0Afzd/cqpzSMpP4kD+AY6XHCfeN56XZ76sy3jD2ZIa9qSVsD+jlP2ZpWSX/TCm4OliJNzPnVBfN0J83Aj0csXX3Rlfd2d83J3wcDHi5mzE3Vn77GwUOBsNOBkNGIXAYED73PqBAIMAgxAYDQIh2r3fyboEhyQlVGTBuX2QtU9LCxacgMaqH65xcqPJN5KTXoP4zsWJQ4YmzrTUkNtSQ/vI4GF0JdQtgABXPwa5+jHIzR8/10F4u/rg5eqDl6sf7s5euDu74+HkgavRFVejK85GZ1yNrjgZnHA2OONkcMIojOqHQx/0NaBPA34tpVzQ+v3TAFLKP7S75m/ANinlv1u/Pw3MllLmdVZubwP6V7te4unU81dEBpokI5pMjGhqYVxDI9NqanBuaTCjNAFeweAV0prvbs1xBw3VfkV1dutx++xdSV0JG89uZG36Wg4XHb7ofSeDE2MDx3JpxKXcPfJu3Jws8wyKqxtILaz+/iOvoo78ygYKKuoprm6guReHVJvLaBAItIkhAi3gaz8MtM9CaBuTtb0mhMDY+rp2nxaMhPhhcklbOee931ZhF7Grq7Cmd9ATgJAmQmQhg1tyCJUFhJnyCTEV4ifL8ZMV+Jkq8KKGBgGZzk6cdXamwGgkz8mJfCcjpUYDZQYjZUYDlQYDspdtNEqJUWoLYQzw/dcCMLT+1bf9Ptj2WWv/D1+3/3zhn7Or7zsjrJh9nuo2jv+9q3cru7sK6OaMhkQAWe2+z0brhXd3TQRwXkAXQiwBlgBERfXumLTEkIm8mr2fQOFKoMGVAKMLXkZXEEatV21wAicXLa/t5KKdseniCc6eWk7bzQ/cB4G7H3gEnjdwNBAEuAdw2/DbuG34bVQ0VFBUW0RJfQkldSUEuAcwNmisVQ7HDvRyJdDLlalxFw8QSympa2qhsq6Zirom6ppaqGtsoa6pmcZmE40tkqZmE80mEy0maJESk0likhIpOe+zqfVzS+v7JpOkpfV9Ca2fte9bTO3v/+G6tvsvvKetCytb29z2XttrbX+WznQZP3QOLvK8An0oJ4Fy4FQH1wppwtVUi7upBndTDS6yATdTHcNM9TjJRpybG3FqasIoG2mW9TRTT5Oop4VGmmmmmSaaaaFZtGCihWZaaBESEy20YMKExASYhAmJRPtv22eNSbR+LSVStP4ddfpoZAdfdXZtZ8+nt8y884KfKoPcgntdY1fMiWYd/YC78E9hzjVIKVcAK0DroZtR90VCE+YTmmDHuyL2I76uvvi6+pJA1/l0axNC4OHihIeLE6G+jvdbkqJYijmjXNlA+6kRkcCFa8vNuUZRFEWxIHMC+gFgiBAiVgjhAtwKfHHBNV8AdwvNVKCiq/y5oiiKor9uUy5SymYhxCPABrRpiyullCeEEEtb318OrEWb4ZKKNm3xHss1WVEURemIWSOCUsq1aEG7/WvL230tgYf1bZqiKIrSE2qzDUVRFAehArqiKIqDUAFdURTFQaiAriiK4iBsttuiEKII6PgsLcsJBNSx6D9Qz+Ni6pmcTz2Pi9n6mURLKYM6esNmAd0WhBBJne2BMBCp53Ex9UzOp57Hxez5maiUi6IoioNQAV1RFMVBDLSAvsLWDbAz6nlcTD2T86nncTG7fSYDKoeuKIriyAZaD11RFMVhqYCuKIriIBwuoAshjEKI74QQX7V+7y+E2CSESGn9PKjdtU8LIVKFEKeFEAts12rLEUJkCiGOCSEOCyGSWl8bsM9ECOEnhPiPEOKUECJZCDFtgD+PYa3/Nto+KoUQ/zPAn8njQogTQojjQoh/CyHc+s3zkK3HbTnKB/Az4CPgq9bv/wg81fr1U8DLrV+PBI4ArkAskAYYbd1+CzyPTCDwgtcG7DMB/gXc3/q1C+A3kJ/HBc/GCOQD0QP1maAdnZkBuLd+/ynwk/7yPByqhy6EiASuAt5p9/K1aP8T0/r5unavfyylbJBSZqDt5T7ZSk21tQH5TIQQPsBM4B8AUspGKWU5A/R5dOByIE1KeZaB/UycAHchhBPggXb6Wr94Hg4V0IHXgF8CpnavhcjW05NaP7edztrZwdaORgIbhRAHWw/phoH7TOKAIuDd1rTcO0IITwbu87jQrcC/W78ekM9ESpkDvAqcQzvkvkJKuZF+8jwcJqALIa4GCqWUB829pYPXHHEO56VSyonAIuBhIcTMLq519GfiBEwE3pZSTgBq0H597oyjP4/vtR4vuRj4rLtLO3jNYZ5Ja278WrT0STjgKYS4s6tbOnjNZs/DYQI6cCmwWAiRCXwMzBVCfAAUCCHCAFo/F7ZePyAOtpZS5rZ+LgT+i/br4EB9JtlAtpRyX+v3/0EL8AP1ebS3CDgkpSxo/X6gPpMrgAwpZZGUsglYBUynnzwPhwnoUsqnpZSRUsoYtF8dv5FS3ol2gPWPWy/7MbCm9esvgFuFEK5CiFhgCLDfys22KCGEpxDCu+1rYD5wnAH6TKSU+UCWEGJY60uXAycZoM/jArfxQ7oFBu4zOQdMFUJ4CCEE2r+RZPrL87D1qLKFRqpn88MslwBgC5DS+tm/3XXPoo1KnwYW2brdFngOcWgj8EeAE8Cz6pkwHkgCjgKrgUED+Xm0/hk9gBLAt91rA/aZAL8BTqF1ft5Hm8HSL56HWvqvKIriIBwm5aIoijLQqYCuKIriIFRAVxRFcRAqoCuKojgIFdAVRVEchAroiqIoDkIFdEVRFAfx/4di14xNteT9AAAAAElFTkSuQmCC\n",
"text/plain": [
"